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I. INTRODUCTION 

The basic ideas of quantum computation were formulated more than 20 years ago 0, Bl- ^^ the 1990's sev- 
eral algorithms have been discovered [21, U, Mi IS l3i IS| that run much faster on a quantum computer than 
on a conventional computer, promising the solution of some problems that are considered to be intractable 
for conventional computers. These discoveries have fueled many new developments, both theoretically and ex- 
perimentally. Quantum information theory is a new interdisciplinary field of research, building on concepts 
from theoretical computer science and theoretical physics |2|. Qn the experimental front, considerable progress 
has been made to engineer and contr ol q uantum systems that may be used for quantum information process- 
ing '11 111 111 111 El til 111 111 111 liilifl III 111 1110 The feasibility of executing 
short quantum algorithms on quantum systems with a few (< 8 bits) has been demonstrated [l4, 15, 16, 17, 29, 30, Slj. 
The technological challenges to build a quantum computer that can perform calculations that would exhaust the re- 
sources of a conventional computer seem tremendous, and there is no indication that such a machine will become 
available in the next few years. 

Any physically realizable quantum computer is a complicated many-body system that interacts with its environment. 
In quantum statistical mechanics and quantum chemistry, it is well known that simulating an interacting quantum 
many-body system becomes exponentially more difficult as the size of the system grows. Actually, it is this observation 
that lead Feynman to the question what kind of computer would be needed to overcome this exponential increase |2] ■ 
Formally, the answer to this question is "a quantum computer." It seems that only a quantum computer can efficiently 
simulate itself |2.l3i^. 

C omp uter simulation has long been accepted as the third methodology in many branches of science and engineer- 
ing |34| . Conventional computers can be used to simulate quantum computers that are relatively small (such as 24 
qubits) but are significantly larger than the experimental machines that have been built. Therefore, it is striking that 
theoretical ideas about quantum computation are seldom confronted with numerical experiments that can be carried 
out on present-day (super) computers. Conventional computers can simulate the abstract model of an ideal quantum 
computer on which most theoretical work is based, and, most important, they can also simulate the physical behavior 
of quantum computer hardware [23 ■ 

There are a number of excellent reviews IHIlllslllllsllillilliil and books 0, H 13 lill that cover the 
theoretical and/or experimental aspects of quantum computations, but there is no review that discusses methods to 
simulate ideal and physical models of quantum computers. The purpose of this review is to fill this void. It contains 
enough information and detail to allow advanced students to develop their own quantum computer simulator. It 
also contains an up-to-date account of existing simulation software. Most of the simulation methods covered in this 
review have a much broader scope than the application to quantum computation might suggest. As a matter of fact, 
much of the impetus to develop these methods stems from problems encountered in the study of quantum dynamics 
of nanomagnets. 

This chapter is organized as follows. Section II gives a brief account of the basics of quantum computation. Section 
III discusses the implemention of simulation models of quantum computers on a conventional computer and also gives 
a survey of the numerical algorithms that are used to simulate these models. A short review of quantum algorithms is 
given in Section IV. Sections V and VI deal with the simulation of ideal and realistic models of quantum computers, 



respectively. In Section VII we give an overview of existing quantum computer simulator software and also provide 
examples of how a simulator is used to perform quantum computation. A brief summary and outlook are provided in 
Section VIII. 

II. THE IDEAL QUANTUM COMPUTER 

In contrast to a digital computer, in which the state of an elementary storage unit is specified by the logical values 
and 1 (= one bit), the state of an elementary storage unit of a quantum computer, the quantum bit or qubit, is 
described by a two-dimensional vector of Euclidean length one. Denoting two orthogonal basis vectors of the two- 
dimensional vector space by |0) and |1), the state 1$) of the qubit can be written as a linear superposition of the 
basis states |0) and |1): 

|$)=ao|0)+ai|l), (1) 

where ag and ai are complex numbers such that |aoP + Wi\'^ = 1- The appearance of complex numbers suggests that 
one qubit can contain an infinite amount of information. However, the principles of quantum mechanics, do not allow 
retrieving all this information ,46- .47- .48.] . The result of inquiring about the state of the qubit, that is a measurement, 
yields either the result or 1. The frequency of obtaining (1) can be estimated by repeated measurement of the 
same state of the qubits and is given by |aop (|aip). 

Just as a conventional computer with only one bit is fairly useless, a quantum computer should also contain more 
than one qubit. Disregarding the enormous technological hurdles that have to be taken to put more than, say 7, 
qubits together, on present-day personal computers we can readily simulate quantum computers with L = 20 qubits. 
The internal state of a quantum computer with L qubits is represented by a unit vector in a Z) = 2^ dimensional 
space (see Section Hi Bp . 

The ideal quantum computer differs from a conventional (probabilistic) computer in the sense that the state of the 
qubits changes according to the rules of quantum mechanics, that is through rotations of the vectors representing 
the state of the qubits ^ . From elementary linear algebra, we know that a rotation of a vector corresponds to the 
multiplication of the vector by a unitary matrix. Thus, the internal state of the quantum computer evolves in time 
according to a sequence of unitary transformations. Any such sequence is called a quantum algorithm. Of course, 
not every sequence corresponds to a meaningful computation. If a quantum algorithm cannot exploit the fact that 
the intermediate state of the quantum computer is described by a linear superposition of basis vectors, it will not be 
faster than its classical counterpart. 

The unitary time evolution of the internal state of the quantum computer is interrupted at the point where we inquire 
about the value of the qubits, that is as soon as we perform a measurement on the qubits. If we perform the readout 
operation on a qubit, we get a definite answer, either or 1, and the information encoded in the superposition is lost. 
The process of measurement cannot be described by a unitary transformation |48| . Therefore, we do not consider 
it to be part of the quantum algorithm. However, to estimate the efficiency of a quantum computation, both the 
operation count of the quantum algorithm and the cost of measuring the outcome of the calculation have to be taken 
into consideration 49]. 

In the quantum computation literature the convention, is to count each application of a unitary transformation as 
one operation on a quantum computer 9] . As the unitary transformation may change all amplitudes simultaneously, a 
quantum computer is a massively parallel machine. To simulate these unitary operations on a conventional computer, 
we actually perform matrix- vector multiplications that requires 0{2^^) (in the worst case) arithmetic operations. 
This may cause some confusion at some points, so the reader should try to keep this in mind. 

To summarize: At any point in time, the internal state of a quantum computer with L qubits is described by a unit 
vector in a Z) = 2^ dimensional space. A quantum algorithm that executes on the quantum computer changes the 
internal state by rotating the vector. Therefore, each step in the quantum algorithm corresponds to the multiplication 
of the vector by a unitary matrix. Readout of the result of a quantum algorithm is destructive in the sense that it 
destroys the information contained in the superposition. 

In the following subsections, we discuss the basic ingredients of quantum computation in more detail. 

A. Spin-1/2 Algebra 

In the simplest form, a qubit is a two-state quantum system, of which there are many examples in the quantum 
world 46, 47, 48] . Measurements of a component of the spin (= internal angular momentum) of particles such as 



electrons, protons, and neutrons along any direction yield either h/2 or ~h/2. The convention is to call the state 
of the spin that corresponds to the outcome ?i/2 "spin up" (| t)) and the other state "spin down" (| |)). The fact 
that the internal angular momentum takes only two values implies that it corresponds to a total angular momentum 
quantum number of 1/2, hence the name spin 1/2 (5 ~ 1/2) particle. 

In principle, any physical object that carries a, S — 1/2 degree of freedom can be used as a physical realization of 
a qubit. In general, the wavefunction that describes the spin of these objects can be written as a linear combination 
of the spin-up and spin-down states [4fil WA li^ : 

|<i>)=ao|T>+ai|i>, (2) 

where oq and ai are complex numbers. It is convenient to normalize the length of the vector |$) to one. Then 
|aop + |ai|2 = l. 

The three components of the spin-1/2 operator S acting on the Hilbert space spanned by the states | t) and | J,) 
are defined by [lllillil 

--K-)' "< ^o)- ''-\{ o-?)' <^' 

in units such that h — 1. According to quantum theory, a measurement of a physical quantity A = A^ of a quantum 
system in a state |<I>) yields a real number, the expectation value |4g.l47ll4q| 

(A) = ($|A|<J>)/($|$). (4) 

From Eqs. (0 and JSJ it follows that 

and we conclude that the representations Eqs. |(5J) and ©, have been chosen such that | t) and | J,) are eigenstates of 
5^ with eigenvalues +1/2 and —1/2, respectively. This is the convention adopted in the physics literature. 
In the context of quantum computation, it is convenient to define [9| 

|0) = |T)=(J) , |i) = U) = (^5). (6) 

Then we have 

|$)=ao|0)+ai|l), (7) 

and the expectation values of the three components of the qubits are defined as 

, , ^ 1 - apg^ - a^ai ^ 1 - zapal + ia^ai ^ 1 - npag + a^Qi 

such that < {Q") < 1 for a = x, y, z. In the following, when we say that the state of a qubit is 0(1), we actually 
mean that (Q^) = 0(1) and (Q^) = {Qv) = 1/2. 

A change of the state of a qubit corresponds to a rotation of the spin. In general, a rotation of the spin by an angle 
(j) about an axis /3 can be written as (?i = 1) 

5" (0, /?) = e*"^-^" S"e-"l"^^ = S" cos cj) + e^f^^S^ sin <j>, (9) 

where use has been made of the commutation rules of the components of the angular momentum operator S |46l . l47ll48| : 

[S'',SP]=ieo^p^S\ (10) 



where tap-j is the totally asymmetric unit tensor {cxyz = ^yzx = ^zxy — 1, ea/37 — —^pa-y = —^-ypa = —^a-yp, 
^aa-y = 0) and the summation convention is used. In addition to Eq. ll()|l the S = 1/2 operators ^ also satisfy the 
relations ^ M, M 

gxgx ^ gygy ^ gzgz ^ 1 (^^^ 

7 11 7 7 ? 

ax qy _ QZ OX oz oy oy qz ^ _ ox oy ox oz oz ox _ oy qz oy __ _ Q^ /"l 9~1 

~ 2 ' ~2' ~2' ~2' ~2' ~2' 

which are often very useful to simplify products of 5 = 1/2 matrices. From Eq. (|12|) it follows that 25*^, 25^, and 
25^ are unitary operators. Hence each of them represents a genuine quantum computer operation. 

Rotations of the spin by 7r/2 about the cc-axis and y-axis are often used as elementary quantum computer operations. 
In matrix notation, they are given by 

x.e'-^'''^±(]l). y..-/=.J=(_; ;), ,,3, 

respectively. X and Y represent operations on single qubits. The matrix expressions for the inverse of the rotations 
X and Y, denoted by X and Y, respectively, are obtained by taking the Hermitian conjugates of the matrices in H13|l . 
With our convention, (0|yiS'^F|0) = —1/2 so that a positive angle corresponds to a rotation in the clockwise 
direction. In general, a rotation of the state about a vector v corresponds to the matrix 

jvs ^ ^ 2iv • S , V ,_, ^, 

e*"^ '^ = 1 cos - + sm-, (14) 

2 w 2 



Vy + v1 is the length of the vector v. 



B. State Representation 

A quantum computer must have more than one qubit. Let us denote by L the number of S* = 1/2 objects of 
the corresponding physical system. In general, the quantum state of a system of L spins is defined by the linear 
superposition of the direct-product states of the L single-spin states: 

1$) = a(tT . . . T)l TT . • • T> + aUT . . . T)UT • . • T> + • • • + a(Ti . . . i)| Ti ■ • • i> + «ai • • • Dili ■ ■ • i>- (15) 

As before, the coefficients (amplitudes) a{]] ... f), ... , a{[[ ■ ■ ■ [) are complex numbers, and it is convenient to nor- 
malize the vector |(/)) by the rescaling 

^ Ha,...aL)? = l, (16) 

0-i...CTi=t,i 

such that (010) — 1. 

In physics, it is customary to let the first label correspond to the first spin, the second to the second spin, and so 
on. Thus, the second term in Eq. H15(l is the contribution of the state with spin 1 down and all other spins up. In 
computer science, it is more natural (and logically equivalent) to think of qubit (spin) 1 as the least significant bit 
of the integer index that runs from zero (all spins up) to 2^ — 1 (all spins down). This is also the convention that is 
adopted in the literature on quantum computation: spin up (down) corresponds to a qubit in state (1). Thus, in 
quantum-computation notation Eq. (|15|l reads 

1$) = a(0 . . . 00)|0 . . . 00) + a(0 . . . 01)|0 . . . 01) + . . . + a(l . . . 10)|1 . . . 10) + a(l . . . 11)|1 . . . 11), 

= ao|0) + ai|l) + . . . + a2^_2|2^ - 2) + a2L_i|2^ - 1). (17) 

For later reference, it is useful to write down explicitly some examples of translating from physics to quantum computer 
notation and vice versa: 



(TTT- 


•T) 


TTT-. 


T) 


= a(0.. 


an- 


•T) 


itt-. 


T) 


= a(0.. 


(Tit- 


•T) 


Tit-. 


T) 


= a(0.. 


(Tii- 


■i) 


Tii-. 


i> 


= a(l.. 


(iii- 


•i) 


iii-. 


i> 


= a(l.. 



.000)10... 000) =ao|0), 
.001)10... 001) =ai|l), 
.010)10... 010) =a2|2), 
.110)|1...110) ==a2i_2|2^-2) 
.lll)|l...lll)=a2^_i|2^-l) 



(18) 



Obviously, there is a one-to-one correspondence between the last line of Eq. 117|l and the way the amphtudes a^ 
are stored in the memory of a conventional computer. From representation H17() . it is easy to estimate the amount of 
computer memory that is needed to simulate a quantum spin system of L spins on a conventional digital computer. 
The dimension D of the Hilbert space (that is the number of amplitudes a^) spanned by the L spin- 1/2 states is 
D = 2^. For applications that require highly optimized code, it is often more efficient to store the real and imaginary 
part of the amplitudes a^ in separate arrays. Furthermore, even for simple looking problems, it is advisable to use 
13 - 15 digit floating-point arithmetic (corresponding to 8 bytes for a real number). Thus, to represent a state of the 
quantum system of L qubits in a conventional, digital computer, we need a least 2^+'* bytes. For example, for L = 20 
(L = 24) we need at least 16 (256) Mb of memory to store a single arbitrary state 1$). Not surprisingly, the amount 
of memory that is required to simulate a quantum system with L spins increases exponentially with the number of 
spins L. 

Operations on the jth spin are denoted by simply attaching the label j to the three components of the spin operator. 
For example, the operation that flips the second spin (or qubit) is given by 



m = 



5^1$) 



a(tTT...T)ITiT...T)+a(iTT...T)liiT-..T) 
+a(Tii.-.i)ITTi...i)+a(iii...i)liTi... 

a(0 . . . 000)|0 . . . 010) + a(0 . . . 001)|0 . . . Oil) + 
+a{l . . . 110)|1 . . . 100) + a(l . . . 111)|1 . . . 101). 



i> 



(19) 



It is important to note that although the spin operator S'^ only flips the jth spin, to compute |$') on a conventional 
computer we have to update all the 2^ amplitudes. In general, on a conventional computer each operation on the 
wavefunction involves at least 2^ arithmetic operations (here and in the sequel, if we count operations on a conventional 
computer, we make no distiction between operations such as add, multiply, or get from and put to memory). The fact 
that both the amount of memory and arithmetic operation count increase exponentially with the number of qubits is 
the major bottleneck for simulating quantum systems (including quantum computers) on a conventional computer. 
However, at the time of writing, the number of qubits that can be simulated far exceeds the number of qubits of 
experimentally realizable systems. 



C. Universal Computer 

It is easy to see that any (Boolean) logic circuit of a conventional computer can be built by an appropriate 
combination of NAND-gates 50] . Although in practice it is often convenient to think in terms of logic circuits that 
can perform more complex tasks than a NAND operation, conceptually it is important to know that NAND-gates are 
all that is needed to build a universal computing device. Also, a quantum computer can be constructed from a set of 
basic gates (= unitary transformations), but instead of one gate we need several |3, ISa, ISM^ IsJ, IS^, IS^ ■ The minimal 
set of gates is not unique, but, for ideal quantum computers, it is mainly a matter of taste to choose a particular 
set. However, in the real world, the kind of operations we can perform on a specific physical system is limited, and 
this will bias the choice of the set of basic gates. Nevertheless, as long as the chosen set of basic gates allows us to 
synthesize any unitary tranformation on the state |$) of the quantum computer, this set can be used for universal 
quantum computation J3, S| • 

In analogy with Boolean logic circuits, simplicity is an important guideline to determine the set of elementary 
unitary gates of a quantum computer. Some basic results from linear algebra are very helpful in this respect. It 
is well known that any vector of N elements can be transformed to the vector (1, 0, . . . , 0) by at most N — 1 plane 
rotations ,.55:. .56] . Each plane rotation is represented by a 2 x 2 unitary matrix that operates on the elements (1, j) for 
j — 2, . . . ,N. Application of this procedure to the A'', mutually orthogonal, column vectors of a A^ x A^ unitary matrix 
U immediately leads to the conclusion that U can be written as a product of at most N{N — l)/2 plane rotations. 



Thus, we can synthesize any N x N unitary matrix by multiplying unitary matrices each of which involves only two 
of the N basis states |57l |. 

As the order in which we apply the plane rotations is arbitrary, the decomposition is not unique. Furthermore, 
on a quantum computer iV = 2^ so that for a given quantum algorithm (or equivalently, a given unitary matrix U) 
we should be able to find a decomposition that takes much less than 2^"-'^ (2^ — 1) plane rotations. Otherwise the 
implementation of the quantum algorithm is not efficient, and there is no point of even contemplating the use of a 
quantum computer to solve the problem. Finding the optimal decomposition oi a N x N unitary matrix U in terms 
of plane rotations is a difficult computational problem, in particular because it is known that for some U there is 
no other alternative than to decompose it into N{N — l)/2 plane rotations Q. It has been shown that any of the 
2^-1^2^ — 1) plane rotations can be constructed by a combination of single qubit gates and the so-called CNOT gate 
that operates on two qubits |3, |^. Therefore these single qubit operations and the CNOT gate constitute a set of 
gates that can be used to construct a universal quantum computer. 

Another basic result of linear algebra is that any nonsingular matrix can be written as the matrix exponential of 
a Hermitian matrix |58l |. As a unitary matrix is nonsingular, we can write U = e^**^, where we already anticipated 
that, for the case at hand, the Hermitian matrix H corresponds to the Hamiltonian that models the qubits. The 
next question is then to ask for the class of model Hamiltonians that can be used for universal computation. Again, 
simplicity is an important criterion to determine useful models, but, as the Hamiltonian represents a physical system, 
there may be some additional constraints on the form of the interactions among qubits and the like. We discuss these 
aspects in more detail in Section IVII where we consider models of physically realizable quantum computers. 

It is a remarkable fact that the simplest quantum many-body system, the Ising spin model, can be used for universal 
quantum computation 51]. The Ising model 

L L 

H{t)^-J2Jl^{t)S^S]-J2 E ^"W^"' (20) 

i,j — 1 i—1 a—x,y,z 

is often used as a physical model of an ideal universal quantum computer 9, '51', '59|. The first sum in Eq. (|20() runs 
over all pairs P < {L— l)L/2 of spins that interact with each other. For instance, if the spins interact only with their 
nearest neighbors, we have P = L/2 and Jij{t) — except when i and j refer to spins that are neighbors. In the case 
of dipolar interaction, P = {L — \)L/2. Note that Eq. (|20|) implicitly assumes that we have complete control over 
all interaction parameters Jf j(t) and external fields hf{t). In Sections III El and III Fl we explain how a universal set 
of gates, the single qubit and the CNOT gates, can be implemented on the quantum computer defined by the Ising 
model lP|l . 

The more general form of the Hamiltonian of a physical model of a quantum computer reads 

L L 

^w = -E E jmsTs<^~Y. E ^"w^"- (21) 

i.j^l a—x.y.z i^l a—x,y,z 

The exchange parameters J" (t) determine the strength of the interactions between the a-components of spins i and 

3- 

Hamiltonian (|21|l is sufficiently generic to represent most models for candidates of physical realizations of quantum 
computer hardware. The spin-spin term in Eq. H21|l is sufficiently general to describe the most common types of 
interactions such as Ising, anisotropic Heisenberg, and dipolar coupling between the spins. Furthermore, if we also use 
spin-1/2 degrees of freedom to represent the environment then, on this level of description, the interaction between 
the quantum computer and its environment is included in model (I21|l . too. In other words, the Hamiltonian 1)21(1 is 
sufficiently generic to cover most cases of current interest. 

In the context of quantum computation and experiments in general, the quantum dynamics of model (|21ll is 
controlled and/or probed by static and/or time-dependent external fields hf{t). Depending on the physical system, 
also the exchange parameters J^j{t) can be controlled by external fields. 

D. Time Evolution of Quantum Computers 

A quantum algorithm is a sequence of unitary transformations that change the state vector of the quantum computer. 
Physically, this change corresponds to the time evolution of a quantum system. In general, the quantum system has 
non-negative probabilities po,Pi, ■ ■ ■ ,P2^-i for being in the states |0), . . . , |2^ — 1). Then, the state of the quantum 
system is described by the density matrix jiy, |43| 



Pit) = J2 \')Mi)('l (22) 

where pi(t — 0) = pi and the expectation values of physical quantities are given by {A{t)) — Trp(t)A. The time 
evolution of the density matrix (|22|l trivially follows from the time evolution of each of the pure states \i). Therefore, 
it is sufficient to focus on methods to compute the time evolution of a pure quantum state. In practice, a calculation 
of the time evolution of the density matrix takes 2^ times more CPU time than the same calculation for a single pure 
state. 

The time evolution of a pure state of the quantum system is given by the solution of the time-dependent Schrodinger 
equation li l liT l li ^ 

i^^mt))=H{t)Mt)). (23) 

The solution of Eq. ^ can be written as [ifiHIvL II^ 

mt + r)) = U{t + T, t)|<i>(i)) = exp+ (-1 ^ H{u)du\ |<i>(t)), (24) 

where r denotes the time step. The propagator U{t + T,t) = exp_|_ [—'iL H{u)du\ is a unitary matrix that 
transforms (that is, rotates) the state |$(i)) into the state \^{t + r)). The time-ordered matrix exponential 
exp_|_ I —i j^ H{u)du\ is formally defined by 

/ rt+T \ rt+T rt+T rui 

exp+ I -i / H{u)du ) = 1 + (-«) / dm H{ui) + {-if / dui / du2H{ui)H{u2) 

rt+T I-Ui /•U2 

+ {-if / dui / du2 / du3H{ui)H{u2)H{u3) + .... (25) 

Jt Jt Jt 

For computational purposes, a naive truncation of the series (|25|l would be fairly useless because it would yield a 
nonunitary approximation to U{t + r, t) and this would violate one of the basic axioms of quantum theory. 

In practice, the time dependence of the external fields can always be regarded as piece- wise constant in time. Hence 
we divide the interval [t,t + t] into n smaller intervals [tk,tk+i] where tk ~ t + X]i=o^i' ''"'= denotes the fc-th time 
interval for which H{t) is constant and t = X]?=o '^i (''o = 0)- Then we can write |6(l| | 

U{t + T,t)^ U{tn, i„-i)C/(i„-i, t„_2) . . . U{ti,to), (26) 

where 

Note that we have not made any assumption about the size of the time intervals [ife,tfe-fi] (with respect to the 
relevant energy scales) . These formal manipulations show that a numerical solution of the time-dependent Schrodinger 
equation for the time dependent model H21(l boils down to the calculation of matrix exponentials of the form H27|l . 
For tj^i <t<tj, the Hamiltonian that appears in H27|) does not change with time. Hence we can use the shorthand 
U{t) = e~**^, keeping in mind that in actual applications, both t and H may depend on the particular time interval 
[tj^i^tj]. 

In general, there are two closely related, strategies to construct algorithms to solve equations such as the time- 
dependent Schrodinger equation [61| . The traditional approach is to discretize (with an increasing level of sophistica- 
tion) the derivative with respect to time. A fundamental difficulty with this approach is that it is hard to construct 
algorithms that preserve the essential character of quantum dynamical time evolution, namely the fact that it is a 
unitary transformation of the state of the quantum system. In particular, in the context of quantum computation it 
is more natural to consider algorithms that yield a unitary time evolution by construction. In this chapter we review 



only methods that are based on the second strategy |6lL l62|. namely, to approximate the formally exact solution, 
that is the matrix exponential U{t) — e^**^ by some unitary time evolution matrix that is also unitary (to machine 
precision) U{t). 

If the approximation U{t) is itself an orthogonal transformation, then ||C/(t)|| = 1 where ||X|| denotes the 2-norm 
of a vector or matrix X [sl HI IH ill ■ This implies that ||*(i)|| = ||c7(i)*(0)|| = ||*(0)||, for an arbitrary initial 
condition $(0) and for all times t. Hence the time integration algorithm defined by U{t) is unconditionally stable by 
construction 61, 62]. In other words, the numerical stability of the algorithm does not depend on the time step t that 
is used. In Sections IIIII we review and compare different algorithms to compute U{t) = e~**^. 

E. Single-qubit Operations 

The simplest, single-qubit operation on qubit j changes the phase of the amplitudes, depending on whether qubit 
j is or 1. In terms of spins, the Hamiltonian that performs this operation is 

H = -h'jS'-, (28) 

and is clearly a special case of the ideal quantum computer model H2U|) . Taking j = 1 as an example, the unitary 
transformation 



t/ = e-**-f^ = e**'*i^i, (29) 



changes the state (|15|l into 



1$') = c/|$) ^ e**a(tT . . . T)l TT • • • T) + e-*^a(it . . . T)UT • • • T> + • • • 

+e^^a(Ti . . . i)| Ti • • • i> + e-*Mii • • • i)Ui . . • i), 

= e'^ [a(TT . . . T)l TT • • • T> + e-''^a(iT . . . T)UT ■ • • T> + ■ ■ • 

+ «(Ti • • • i)l Ti . . • i> + e-'^^a(ii . . . i)Ui . . . i)] , 

3*"* [ao|0) + e-2'^ai|l) + . . . + a^L^^l^^^ - 2) + e-2^%._i|2^ - 1)] , (30) 



e 



where the phase shift (j) is given by = i/if /2, and we used the fact that in the (spin up, spin down)-representation, 
Sf is a diagonal matrix with eigenvalues (1/2,-1/2). From Eq. (jHOfl it is clear that by a proper choice of the time t 
and external field hf we can perform any phase shift operation. To obtain the last line of Eq. (|30|l . a global phase 
factor has been extracted from all the amplitudes. This manipulation does not alter the outcome of the calculation 
as it eventually drops out in the calculation of expectation values |4y, l47l l4^ . We will often use this trick to simplify 
the expressions, dropping global phase factors whenever possible. 

The phase shift operation Zj = e^'^^i with the 7r/2 rotations (|13|) allows us to compose any single qubit operation. 
For example, using the algebraic properties of the spin operators, it follows that 

showing that we can readily construct a 7r/2 rotation about the z-axis by combining rotations about the x and y-axis. 
Some simple examples may help to understand the various operations. For a two-qubit quantum computer we have 
\a) = ao\00) + ai|01) + azIlO) -|- aajll) and 
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(32) 



where Xi denotes a rotation of qubit 1 by tt/2 about the x-axis. For example Xi|10) = (|10) + i\ll))/V^ and 
Xi|10) = (|10) — i|ll))/-\/2. Using the same labeling of states as in (pT2)l . we have 



TABLE I: Input and output states and the corresponding expectation values (QijQI) of the qubits for the CNOT operation. 



Input state 



Output state 



100) 
101) 
110) 

111) 



OjOO) 
0111) 

i|io) 

1101) 



'--T^ 



/ 1 1 

10 1 

-10 10 

V -1 1 



(33) 



for example, FallO) = (|00) + |10))/V2 and FallO) = (-|00) + |ll))/\/2. 
For later use we introduce the symbol 



Rj{(j))=e'^'^e-''*'^^ = ( 



1 

e**^ ' ' 



(34) 



to represent the single-qubit phase-shift operation by a phase on qubit j. A symbolic representation of Rj{4) — 
Ti/k) is given in Fig. ^. 

F. Two-qubit Operations 

As explained previously, the single qubit operations and the CNOT gate constitute a universal set of gates 9] . The 
CNOT gate is defined by its action on the computational basis states, as shown in Tabled To simplify the notation, 
we consider only the two relevant qubits and call them qubit 1 and 2. The CNOT gate flips the second qubit if the 
state of the first qubit is |1), that is, the first qubit acts as a control qubit for the second one. As shown in Fig.^, the 
CNOT gate can be symbolically represented by a vertical line connecting a dot (control bit) and a cross (target bit). 
On the ideal quantum computer model H20|l . the CNOT gate can be implemented by a combination of single-qubit 
operations and a controlled phase shift operation. In matrix notation, the CNOT operation (see Table QJ can be 
written as 
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Our goal is to show that up to an irrelevant global phase factor, CNOT — 1^2-^21 (7r)F2 where 



(35) 
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(36) 



It is easy to see that 



hi{^) 



t{-JSlS^-hSl-hS^) 



(37) 



where h — —J/2 and Jt — (j). Eq. I|37(l shows that physically the phase-shift operation H36|l corresponds to the time 
evolution of a quantum spin system described by the Ising model. 



Using the properties of the S = 1/2 matrices, it is easy to show that 



CNOT = Y2 
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y2 = e-"/4r2/2i(vr)r2. 



(38) 



This completes the construction of the set of universal gates that can be implemented on the ideal Ising-model 
quantum computer (|20|l . 

The CNOT operation is a special case of the controlled phase shift operation Rji{(^). The latter appears in 
several of the quantum algorithms that are discussed later and therefore it is appropriate to introduce it here. The 
controlled phase shift operation on qubits 1 and 2 reads 



i?2i(0) = e'^/"/2i(<^) 



10 
10 
10 

e^"^ . 



(39) 



Graphically, the controlled phase shift i?^ (0 — ii/k) is represented by a vertical line connecting a dot (control bit) 
and a box denoting a single qubit phase shift by -n/k (see Fig.^^). 
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FIG. 1: Graphical representation of some of the basic gates used in quantum computation; (a) single qubit phase shift 
Rj{(t) = 7r/fc) (see Eq.OU; (b) CNOT gate (see Eq.OHJ; (c) controlled phase shift by _Rji(<^ = 7r/fc) (see Eq.EHJ; (d) Walsh- 
Hadamard gate (see Eg. 1761 : (e) SWAP gate (see Ea. l78ll : (f) Toffoli gate (see Ea. l83ll . The horizontal lines denote the qubits 
involved in the quantum operations. The dots and crosses denote the control and target bits, respectively. 



III. NUMERICAL METHODS 



We start by discussing aspects that are common to all numerical algorithms discussed in this review such as 
techniques to let (combinations of) spin-1/2 operators act on a state. Then we review four different algorithms to 
compute the time evolution C/(i)|$) — e~'*^|$). We start by a brief exposition of the fuU-diagonalization approach. 
Then we describe the Chebyshev polynomial algorithm, the short-iterative Lanczos procedure, and several algorithms 
based on Suzuki- Trotter product formulas. 

In numerical work where the use of floating-point arithmetic (e.g. 13-15 digit precision) is necessary, there always 
will be errors due to the finite precision, rounding and the like. In this chapter we use the term numerically exact 
to indicate that the results are correct up to an error that vanishes with the number of digits that are used to perform 
the calculation. For example, the value of arccos(— 1) = tt is numerically exact. 



A. General Aspects 

As explained previously, the state of a spin- 1/2 system with P pairs of L interacting spins or, equivalently, the state 
of a L-qubit quantum computer, is represented by a complex- valued vector of length D = 2^ . In quantum mechanics 
matrices operating on a vector of length D have dimension D x D. Matrix elements of the Hamiltonian (|21|l can 
be complex numbers. Hence we need 2^^"'"'* bytes to store the full matrix. For instance, in the case of a 10-qubit 
quantum computer we need 16 Mb (16 kb) of memory to store the Hamiltonian (state) which poses no problem 
for present-day PCs. However, to store the Hamiltonian (state) of the L = 20 system, a computer should have 16 
Tb (16Mb) of memory. This simple example clearly shows that any simulation scheme that requires storage for the 
matrix representing the Hamiltonian will be of (very) limited use. Thus, it is necessary to consider techniques to 
avoid the use of large matrices. Fortunately, for the spin systems considered in this review it is rather straightforward 
to organize the calculations such that we need only storage for 0{PL) matrices of dimension 2x2 and 4x4. 

We first consider the operations |$') = S',"|$). More specifically, let us denote 

1$) = a(TT . . . T)l TT . • • T> + a(iT . . . T)UT • ■ • T> + • • ■ + a(Ti . . . i)| Ti • • • i) + «(ii • • 4)Ui ■ • • i), (40) 

and 

1$') = Sf\^) 

- a'in . . . T)l TT • • . T) + a'(iT . . . T)UT • . • T> + • . • + a'(TT . . . i)| Ti ■ • • i) + a'Ui . . . T)Ui . . . i). (4i) 

From Eq. 0, it follows that Sj is a diagonal matrix in the representation that we use to store the vector |$). Thus 
we obtain |$') by reversing the sign of all coefficients of |$) for which the jth bit of their index is one. In terms of 
amplitudes, we have 

a'(* ... * * ... *) = +— a{* ... * * ... *) 

a'{* . . . * I * . . .*) — — a{* . . .* 1 * . . .*), (42) 

where we use the asterisk to indicate that the bits on the corresponding position are the same. From Eq. Q, it follows 
that Sj interchanges states up and down. It immediately follows that the elements of 1$') are obtained by swapping 
pairs of elements of | $) . We have 

a'{* ... * * ... *) = +— a{* ... * 1 * ... *) 

a'{* . . . * I * . . .*) = +- a{* . . . * * . . .*). (43) 

In Eq. H43|l . the bit strings on the left- and right-hand side are identical except for the jth bit. In the case of 
1$') = S'j'l^), a similar argument shows that 

a'(* ... * * ... *) — — a{* ... * 1 * ... *) 

a'{* . . . * 1 * . . .*) ~ +-a{*...*0*...*). (44) 

Note that all these operations can be done in place, that is, without using another vector of length 2^. The three 
operations discussed here are sufficient to construct any algorithm to simulate the quantum dynamics of model H21|l . 
However, for reasons of computational efhciency, it is advantageous to extend the repertoire of elementary operations 
a little. 

The two-spin operation |$") = SfS^l^) can be implemented in at least three different ways. In the discussion 
that follows, we exclude the trivial case where j = k. One obvious method would be to write it as two single-spin 
operations of the kind discussed previously: |$') — S^\^) followed by 1$") — S"\^'). For a single pair {j, k), this is an 
efficient method. However, if there are many pairs, as in the case of Hamiltonian (|21() . we have either to recalculate 



S'^|<I>) several times, or we have to allocate extra storage to hold S'^|<i>) for fc = 1, . . . , L. In both cases, we use extra, 
costly resources. 

In the second approach, we work out analytically how the amplitudes change if we apply S'^S^ to |<I>). li a — z, 
we have to reverse the sign of the amplitude only if the jth and fcth bit of the amplitude index are different. We have 

a"{* ...*0*...*0*...*) = +— a{* ...*0*...*0*...*) 
a"{* ...*1*...*0*...*) = — a{* ...*1*...*0*...*) 
a"{* ...*0*. ..*!*...*) — — a{* ...*0*. ..*!*...*) 

a"{* ...*!*... *!*...*) = +— a{* ...*!*...*!*...*). (45) 

li a ~ X, we interchange quadruples of amplitudes according to the following, rather obvious, rules: 

a"{* ...*0*...*0*...*) — — a(* ...*1*. ..*!*...*) 

a"{* ...*1*...*0*...*) = — a{* ...*0*. ..*!*...*) 

a"{* ...*0*. ..*!*...*) — — a(* ...*1*...*0*...*) 

a"{* ...*!*...*!*...*) — — a(* ... * * ... * * ... *). (46) 

We leave the case a = y as an exercise for the reader. Clearly, this approach does not require additional storage or 
calculations to compute intermediate results. 

The third method relies on the observation that for a — x,y and j ^ k 

SfS^m = R'^RtS]St{R'^)HRt)^<i>), (47) 

where R" is a rotation that transforms S^ into S" and that the calculation of S^S^\^) can be done very efficiently. 
From Eq. @ it follows that YjS^Yj = 5| and XjS^Xj = S^ . Hence it follows that i?J = Yj and R^. = Xj. From 
Eq. (|47|l . it is clear that we have to determine the rules only to compute 1$') — (i?^)'f|<I>) for a — x,y. These rules 
follow directly from the matrix representations H13|l of Xj and Yj. For instance, in the case of |$') = (i?p^|$), we 
have 

a'{* ... * * ... *) = H — p[a(* ... * *...*) — a{* ...*!*...*)] 
v2 

a'{* ...*!*...*) = +—=[a{* ... * *...*) + a{* ...*!*... *)], (48) 

v2 

and similar expressions hold for the other cases. 

All the single- and two-spin operations discussed here can be carried out in 0(2^) floating-point operations. The 
full diagonalization method, the Chebyshev and short-iterative Lanczos algorithm to be discussed later require a 
procedure that computes |$') = H\^). Given the parameters of H [see Eq. (|21(l ]. it is straightforward to combine the 
single- and two-spin operations that were described to perform the operation 1$') = ^I'l')- 

The single- and two-qubit operations described before suffice to implement any unitary transformation on the state 
vector. As a matter of fact, the simple rules described are all that we need to simulate ideal quantum computers on 
a conventional computer. In Section we discuss this topic in more detail. For now, we continue with reviewing 
numerical techniques for computing the time evolution e~**^ of a physical system described by a Hamiltonian H. 

B. Full Diagonalization Approach 

As H is a D X D Hcrmitian matrix, it has a complete set of eigenvectors and real- valued eigenvalues |55l ISm . IsM \6^ . 
Let us denote the diagonal matrix of all these eigenvalues by A and the unitary matrix of all these eigenvectors by V. 



Then we have V^ HV = A and U{t) — e^**'^ — Ve~'^^^V'^ . In other words, once we know A and V , U{t) is obtained 
by simple matrix multipUcation. Thus, the most straightforward method to compute U{t) — e^'*^ = Ve~^^^V^ 
is to use standard hnear-algebra algorithms to diagonalize the matrix H . An appealing feature of this approach is 
that the most complicated part of the algorithm, the diagonalization of the matrix H , is relegated to well-developed 
linear algebra packages. Memory and CPU time of the standard diagonalization algorithms scale as D^ and D^ , 
respectively .55, 56, 63]. 

The matrix elements of H can be computed by repeated use of the |<i>') — H\^) procedure. If we set |<I>) — 
(1, 0, ... , 0)"^, then |$') = H\^) gives us the first column of the matrix H . The same calculation for |<I>) = (0, 1, . . . , 0)^ 
yields the second column, and so on. 

The application of the fuU-diagonalization approach is limited by the memory and CPU resources it requires. The 
former is often the most limiting factor (see the prior discussion) , in spite of the fact that full diagonalization takes 
0{2^^) floating-point operations. In practice, solving the time-dependent Schrodinger equation for a system of 12- 
13 spins for many pairs {H ,t) (which is necessary if there are time-dependent fields) is many orders of magnitude 
more expensive in terms of computational resources than solving the same problem by the methods discussed later. 
Nevertheless, any toolbox for simulating quantum spin dynamics should include code that is based on the full diago- 
nalization approach. For all but the most trivial problems, it is an essential tool to validate the correctness of other 
algorithms that solve the time-dependent Schrodinger equation. 

C. Chebyshev Polynomial Algorithm 

The basic idea of this approach is to make use of a numerically exact polynomial approximation to the matrix 
exponential U{t) — e~'*^ J£J, l^a, Ifia H3, Isl, IS^, llS ■ The first step in the algorithm is to "normahze" the matrix H 
such that its (unknown) eigenvalues lie in the interval [—1,1]. In general, the eigenvalues Ej of the Hermitian matrix 
H are real numbers in the interval [— ]]_ffj|, \\H\W where \\H\\ = maxj \Ej\ is the 2-norm of the matrix H 55^^]. Unless 
we already solved the eigenvalue problem oi H, we usually do not know the exact value of \\H\\. For the Hamiltonian 
()21|l . it is easy to compute an upper bound to ||_ff ]] by repeated use of the triangle inequality 

l|A: + r|| < l|xi| + ||yi|, (49) 

and the elementary bounds WSfW = 1/2 and US'" 5*^11 = 1/4. For fixed t, we can write H — H{t), and we find |23 

\\H\\<\\H\u^lj2 E \J"J + IJ2 E i^"I- (50) 

i.j — l a—x.y,z i—1 a—x.y,z 

Maximum efficiency of the Chebyshev polynomial algorithm is obtained if we know the exact value \\H\\ (see later). 
Thus, if we have more specific information about the model parameters J"- and hf, it should be used to improve the 

bound on \\H\\. By construction, the eigenvalues of H = H/\\H\\i, all lie in the interval [—1, 1]. 
Expanding the initial state <I'(0) in the (unknown) eigenvectors \Ej) of H, we have 

lci>(t)) = e-*«l$(0)) = e-'^l$(0)) =Y,e-''^^E,){E,mO)), (51) 

j 

where t = iH^Hh and the Ej — Ej/\\H\\i, denote the (unknown) eigenvalues of H (in practice, there is no need to 
know the eigenvalues and eigenvectors of H explicitly). Next, recall that for \x\ < 1 we have |7]| 

oo 

e-'^^- Jo(z) + 2^(-z)'=Jfe(2)Tfe(x), (52) 

fe=i 

where Jfc(z) is the Bessel function of integer order k, and Tk{x) = cos[k arccos(a;)] is the Chebyshev polynomial of the 
first kind 0. As — 1 < i?j < 1, we can use Eq. l(5^ to write Eq. l(5T|l as 



mt)) = 



Jo(-f)]l + 2^Jfe(-t)ffc(ij) 



fe=i 



1$(0)), (53) 
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FIG. 2: Dependence of the Bessel function Jk{z) on the order k; sohd hne: z 
z = 2000. 
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where Tk{H) = i^TkiH) = i^Tlj \^3)T^k{Ej){Ej\ is a matrix- valued, modified Chebyshev polynomial that is defined 
by the recursion relations 



and 



ro(ff)|*)-|*), T^{H)\^) = iH\^) , 



Tk+i{H)\^) = 2iHTk{H)\^) +Tk-i{H)\^) 



(54) 



(55) 



for fc > 1. _ _ 

For practical purposes the sum in Ea. H53() has to be truncated. We compute the approximation |$(i)) ~ f7(t)|$(0)) 
by keeping the first K + 1 contributions: 



l$(i)> = 



K 



M-t)t + 2Y,M~t)Tk{H) 



k=l 



|$(0)). 



(56) 



As I Ja;(— z)| < 1 and ||rfe(iJ)|j — \\Tk{H)\\ < 1, all contributions in Eq. (|56|l are of the same order of magnitude. In 
contrast, the truncated Taylor series e~**^ w X]fe=i(~*^^)'^/^' i^ ^ ^u™^ ^^ many small and large numbers, and a 
numerical calculation of the sum is known to suffer from severe instabilities 72| . 

Using the downward recursion relation of the Bessel functions, it is straightforward to calculate the first K Bessel 
functions to machine precision in 0{K) arithmetic operations |7ll l73|. In practice, a calculation of the first 20,000 
Bessel functions takes less than 1 second on a Pentium III 600 MHz mobile processor, using 14 - 15 digit arithmetic. 
To get a feeling for the value of K at which the Chebyshev polynomial expansion might be truncated, it is instructive 
to plot Jfc(z) as a function of fc. From Fig[21it is clear that | Jfc(z)| vanishes (very) rapidly if fc becomes larger than z. 
In fact, I Jfc(-z)| < z^ /2^k\ for real z hence | Jfc(2:)| vanishes exponentially fast as fc increases [ll|. Thus we can fix the 
number K by requiring that | Jfc(^:)| > k for all k < K. Here k is a small number that determines the accuracy of the 
approximation. From Fig|2 it follows that K ^ z and that \Jk{z = 2000)| < 10^^° for aU fc > z + 100. 

Once we have found the smallest K such that \Jk{—t)\ > n for all k < K, there is no point of taking more 
than K terms in the expansion. Since |jT'fc(i?)|| < 1, Fig|21 suggests that including such contributions would only 
add to the noise, and numerical tests show that this is indeed the case. However, taking less than K terms has 
considerable negative impact on the accuracy of the results. Hence in practice the choice of K is rather limited (such 
as K fn z -\- 200 if z = 2000 and machine precision is required). Most important, for fixed k, K increases linearly with 
t||i7||b. Therefore, if we can replace the bound ||-ff||b by a sharper one, the value of K can be reduced accordingly 
and the calculation will take less CPU time. 



In a strict sense, the truncated Chebyshev polynomial in Eq. (|56|l is not a unitary matrix. However, the error ||$(t) — 
^{t)\\ between the exact state $(t) and the truncated Chebyshev polynomial result <i>(t) vanishes (exponentially) fast 
as K increases. Therefore, for practical purposes the truncated Chebyshev polynomial can be viewed as an extremely 
stable time- integration algorithm because it yields an approximation to the exact time evolution operator U{t) = e^'*^ 
that is numerically exact. Under these conditions, the Chebyshev algorithm may safely be used repeatedly to perform 
multiple time steps with a (very) large fixed time step 70] . 

According to Ea. H56|l . performing one time step amounts to repeatedly using recursion relation H55|l to obtain 
Tk{H)^{Q) for k = 2,...,K, then multiply the elements of this vector by Jk{—t) and add all contributions. This 
procedure requires storage for two vectors of the same length as $ and a procedure that returns H\^). Both memory 
and CPU time of this method are 0{D = 2^). CPU time increases linearly with the time step t. 

D. Short-Iterative Lanczos Algorithm 

The short-iterative Lanczos algorithm |66l ItJL rZa . VWL l77l| belongs to the class of reduced-order methods that are 
based on the approximation of a matrix by its projection onto a (much) smaller subspace. The short-iterative Lanczos 
algorithm is based on the approximation 

e-**^|*) w UN{t) = e-**^"^^"|*), (57) 

where Pjsi is the projector on the iV-dimensional subspace spanned by the vectors {^, i/vj/^ . . . ^ i/^^i\]/}. As 
^-iiPnHPn jg unitary the method is unconditionally stable. 

There are two major steps in the short-iterative Lanczos algorithm. First we calculate PnHPn'^ by generating the 
orthogonal Lanczos vectors in the usual manner 5^ l56] : 

l*o) = I*), 

l*o) = \K)/^J{%\%), 

|*'i) = i/|*o)-|*o)(*o|i^|*o), 



1*1) = \^>[)/^{^i>[\^>[), (58) 

and 

|*'fe) = i7|*fc_i)-|*fc_i)(*fc_i|ff|*fc_i)-|*fc_2)(*fc-2|^|*fc-2), 



\^k) = \^',) /,/{%]%), (59) 



^N 



for 1 < fc < iV. By construction, we have {^k\'^k') = Sk.k', and the projection with Pn = J2k=i l^'s)(^fc| yields a 
N X N tri-diagonal matrix Pi\jHPi\j. The second step in the short-iterative Lanczos algorithm is to diagonalize this 
tridiagonal matrix and use the resulting eigenvalues and eigenvectors to compute e"**-^™^-^" |*). The latter is done 
in exactly the same manner as in the case of the full diagonalization method. 

The short-iterative Lanczos algorithm requires storage for all the eigenvectors and/or the N Lanczos vectors |*fe) 
of the N X N matrix PmHPm- Thus, for this method to compete with the full diagonalization method, we require 
N <S^ D — 2^. The accuracy of this algorithm depends both on the order N and the state |*) |6ftl74l l75|. With infinite- 
precision arithmetic, e~'*^|\l/) = Muin^oo e~**^"^-^™ |\1/), but in practice, the loss of orthogonality during the Lanczos 
procedure limits the order A^ and the time step t that can be used without introducing spurious eigenvalues |55l ISq . 
The low-order short-iterative Lanczos algorithm may work well if |*) contains contributions from eigenstates of H that 
are close in energy. However, if |*) contains contributions from many eigenstates of H with very different energies, 
it is unlikely that all these eigenvalues will be present in PmHPm, and the approximation U]s[{t) may be rather poor 
unless A^ is sufficiently large. Memory and CPU time (for one time step) of the short-iterative Lanczos algorithm 
scale as D and N^D, respectively. In general, A^ increases with f in a nontrivial, problem-dependent manner that 
seems difficult to determine in advance. 



E. Suzuki- Trotter Product-Formula Algorithms 

A systematic approach to construct unitary approximations to unitary matrix exponentials is to make use of the 
Lie- Trotter-Suzuki product-formula |78ll79l | 

I K \ ™ 

[/(t) = e-*'^ = e-'*(^i+-+^-) = lim FT e-^*^'./™ , (60) 

\fc=l / 

and generalizations thereof [Sfl Isil l82 | . Expression H6()(l suggests that 

t7i(t) = e-**-f^i...e-^*^^, (61) 

might be a good approximation to U(t) if i is sufficiently small. From Eq. H61|l . it follows that the Taylor series of 
U(t) and U\(i) are identical up to first order in t. We call U\(f) a first-order approximation to C/(t). If all the Hi 
in Eq. I|61() are Hermitian, then t/i(t) is unitary by construction, and a numerical algorithm based on (|61|l will be 
unconditionally stable. For unitary matrices Uif) and Ui{t), it can be shown that |62l | 



f/W-C/i(i)||<yEll[^-^^]ll' (62) 



i<j 



suggesting that Ui (t) may be a good aproximation if we advance the state of the quantum system by small time steps 
t such that t\\H\\ <C 1. Note that this is the situation of interest if there are external time-dependent fields. 

The Suzuki- Trotter product-formula approach provides a simple, systematic framework to improve the accuracy of 
the approximation to U{t) with marginal programming effort and without changing its fundamental properties. For 
example, the matrix 

U2{t) = Ul{-t/2)Ui{t/2) = e-''""''^ . . .e-**-f^i/2e-**-f^i/2 . . .e-''""''^ , (63) 

is a second-order approximation to U{t) j^OjISHSj- ^^ Ui{t) is unitary, so is U2{t). For unitary [/2(i), we have [S2| 

\\U{t)-U2{t)\\<C2t\ (64) 

where C2 is a positive constant |6^ . 

Higher-order approximations based on U2{t) [or Ui{t)\ can be constructed by using Suzuki's fractal decomposi- 
tion 80, 82] . A particularly useful fourth-order approximation is given by 8ii 82j 

U^it) = U2{at)U2{at)U2{{l - Aa)t)U2{at)U2{at) , (65) 

where a = 1/(4 — 4^^/"^). As before, if U2{i) is unitary, so is UA{t) and we have 

\\U{t)~Ui{t)\\<Cit\ (66) 

where C4 is a positive constant. The rigorous error bounds (|62|l . (|64|l and H66|) suggest that for sufficiently small t, the 
numerical error ||C/(t)^ — C/„(i)^|l vanishes as t". If this behavior is not observed, there is a fair chance that there is 
at least one mistake in the computer program that implements Un{t). For a fixed accuracy, memory and CPU time 
of an Unit) algorithm scale as 0{D) and C'(t'^^+^/"^I?), respectivel y. T he approximation s l|63ll a nd 116511 have proven 
to be very useful for a wide range of different applications |H1 IH IsT, "s?, 84, "s?, "s?, "s?, 88, 89, "9^, 91, "9?, 93, "9^. 

In practice, an efficient implementation of the first-order scheme is all that is needed to build higher-order algorithms 
(|63ll and H65|) . The crucial step that we have to consider now is how to choose the Hermitian Hi& such that the matrix 
exponentials exp(— iiT^i), ..., exp(—itHK) can be calculated efficiently. 

If there are external time-dependent fields, as in the case of a physical quantum computer, it is expedient to 
decompose the Hamiltonian into single-spin and two-spin contributions. The rational behind this is that in most 
cases of interest (NMR,...) the external fields change on a much shorter time scale than the other parameters in the 
Hamiltonian. As the spin operators with different qubit labels commute, we have 
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(67) 



The jth factor of the product in 167|l rotates spin j about the vector hj = {h^,h^-,h^). Each factor can be worked 
out analyticaUy, yielding 
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(68) 



where hj = \\h.j\\ is the length of the vector hj. From Eq. (|67|l and Eq. (|68() . we conclude that the time evolution due 
to single-spin terms can be calculated by performing 0{LD) operations involving 2x2 matrices. Because all matrix 
elements are non-zero, these single-spin operations are only marginally more complicated than the ones discussed in 
Section ImH 

We now consider two different decompositions of the two-spin terms that can be implemented very efficiently: The 
original pair-product split-up |79l l9fT| in which Hj contains all contributions of a particular pair of spins, and an 
XYZ decomposition in which we break up the Hamiltonian according to the x, y, and z components of the spin 
operators |93|- 

The pair-product decomposition is defined by [23,123 
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where the order of the factors in the product over (j, k) is arbitrary. This freedom may be exploited to increase the 
execution speed by vectorization and/or parallelization techniques. Also in this case the expression of each factor can 
be worked out analytically and reads 
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(70) 



where a = Jlk/^, b = {J^,^ - Jlk)/^, and c = (Jl;^ + ^Ik)/^- F™™ Eq. ^ and Eq. ^, it follows that the time 
evolution due to the spin-spin coupling terms can be calculated by performing 0{PD) operations involving 4x4 
matrices. These 4x4 matrix operations are marginally more complicated than the ones discussed in Section IlII Al 
and therefore we omit further details. 

The XYZ decomposition is defined by [23 
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(71) 



where there is no particular order of the factors in Eq. (|71|l and 



iJ" = - 



3.k=l 



ja qq. QCt 



(72) 



The computational basis states are eigenstates of the 5^ operators. Thus, in this representation e^ is diagonal by 
construction, and it changes the input state by altering the phase of each of the basis vectors. As H^ is a sum of pair 
interactions, it is trivial to implement this operation as a sequence of multiplications by 4 x 4 diagonal matrices. Still 
working in the same representation, an efficient algorithm that implements the action of e^' (e^' ) by using the 
Yj (Xj) rotations can be constructed as follows. Writing Y = Y[i=i ^j ^^ have 



g ^rtH- ^ YYe-'*"^YY = Ycxp[itJ2 Jj.kSjS'k ] Y- (73) 

i,k=l 

From Eq. H73|) . it is clear that the action of e~**'^ can be computed by applying the L single-spin rotations Y , a 
sequence of multiplications by 4 x 4 diagonal matrices containing phase shifts, and the L inverse rotations Y . A 
similar procedure is used to compute the action of e~'*-^ . We have only to replace Y (Jf^) by X {J^ j.)- The total 
arithmetic operation count of the XYZ decomposition (|71|l is 0{LD). If there is enough memory on the conventional 
computer to store all the values of Ylj k=i J^,k'^j'^k ^^i' o/. = x,y,z and Sj = ±1/2, then the XYZ decomposition runs 
in 0{D). 

F. Comments 

All the algorithms just discussed satisfy the basic requirement (unitarity to machine precision) of a valid quantum 
mechanical time evolution. A general analysis of the strong and weak points of these algorithms is rather difficult: 
experience shows that the conclusions may change significantly with the particular application. Therefore, we can 
only consider specific examples and compare algorithms in terms of accuracy, memory, and CPU requirements, and 
that is what we will do in this subsection. 

As an example, we consider a model of two spins (Si, S2) interacting with a "bath" of L — 2 spins (I„) described 
by the Hamiltonian |9fij 

L-2 

H = Jo(Si + 82)^ + Y. J„I„ • (Si + S2). (74) 

n=l 

This model has been used to study decoherence in a two-spin system due to coupling with a bath of spins |2a|- 
The Heisenberg exchange interactions {Jn\ are assumed to be random, uncorrelated, and small compared with Jp 
{\Jn\ ^ \Jo\)- Initially, the two spins Si and S2 are in the state with one spin up and the other spin down. The initial 
state of the spins {In} is assumed to be random. Obviously, Eq. (|74|l is a special case of the generic Hamiltonian H21|l . 

In Fig. 13 we show simulation results of the time evolution of the magnetization {Sf{t)) over a large time span, 
for a spin-bath containing 22 spins (L — 24). On a fine time scale (compared with the scale used in Fig|31l, the first 
and second spin oscillate rapidly with a period of order 1/ Jq (results not shown) 96] . Initially, the amplitude of the 
magnetization oscillations rapidly decays to zero, then increases again, and then decays to zero very slowly '96]. This 
is the generic behavior of the magnetization in this class of models ^9&\ . 

We use the system (|74|) to compare the different time-integration algorithms. In Table ^ we show the error on 
the final state and the CPU time it took to solve the time-dependent Schrodinger equation. The error is defined as 
]]^*(to) — \E';j^(m)]] where * denotes the algorithm (exact diagonalization or Chebyshev polynomial method) that is 
used to generate the reference solution and X is one of the other algorithms. It is clear that the short iterative Lanczos 
method is not competitive for this type of time-dependent Schrodinger problem. The fourth-order pair-approximation 
is close but is still less efficient than the Chebyshev algorithm. The other Suzuki product-formula algorithms are 
clearly not competitive. The reason that the pair-approximation is performing fairly well in this case is related to the 
form of the Hamiltonian H74|) . The present results support our earlier finding 70] that the numerical simulation of 
decoherence in spin systems is most efficiently done in a two-step process: The Chebyshev algorithm can be used to 
make a big leap in time, followed by a Suzuki product-formula algorithm calculation to study the time dependence 
on a more detailed level. 

In Tabic IIIII we present a necessarily rough overall assessment of the virtues and limitations of the algorithms 
discussed in this review. From a general perspective, to increase the confidence in numerical simulation results, it is 
always good to have several different algorithms to perform the same task. 

IV. QUANTUM ALGORITHMS 

A quantum algorithm is a sequence of unitary operations that changes the internal state of the quantum computer J9( . 
Quantum algorithms vary in complexity. In the first subsection, we discuss the elementary gates of a quantum 
computer. Two of these gates - namely the CNOT gate and the controlled phase shift - have already been discussed 
in Section IIIFI These quantum gates are rather simple quantum algorithms that can be used to build quantum 
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FIG. 3: Left: Magnetization (Sf (t)) as a function of time as obtained by numerical simulation of two spins interacting with 
a bath of 22 spins. The parameters of model 1741 are Jo = 8 and the J„ are uncorrelated random numbers that satisfy 
< Jn < 0.4. Right: log KS'i (f))| as a function of time showing exponential decay at long times (note that the intrinsic time 
scale of 5*1 (i) is set by Jo = 8). 



TABLE IL Comparison of algorithms to solve the time-dependent Schrodinger equation for model l|74^ . All calculations use a 
time step r/27r = 0.01. L is the total number of spins, m the number of time steps, and N denotes the number of iterations 
in the short-iterative Lanczos algorithm. An asterisk in a column indicates that the state obtained by that particular method 
is used as reference to compare with states obtained by other algorithms. The calculations for L = 22 where carried out on a 
Cray SVl system. All other calculations were performed on an IBM T40 Notebook with an Intel Centrino (1.3GHz) processor 
and 512 MB memory. 





Exact 


Chebyshev 


Iterative 


Suzuki 


Suzuki 


Suzuki 


Suzuki 




Diagonalization 


Polynomial 


Lanczos 


Pair(2) 


Pah(4) 


XYZ(2) 


XYZ(4) 




L = 10, m = 400 
















Error (N=5) 


* 


0.34E-12 


0.17E-05 


0.23E-03 


0.75E-08 


0.14E+00 


0.53E-04 


CPU time [s] 


72 


0.5 


3.5 


0.8 


3.6 


0.6 


2.6 


L = 12, m = 400 
















Error (N=5) 


- 


* 


0.27E-05 


0.27E-03 


0.80E-08 


0.14E+00 


0.55E-04 


CPU time [s] 


- 


2.1 


17.3 


3.7 


18.0 


2.5 


12.1 


L = 12, m = 400 
















Error (N=10) 


- 


* 


0.81E-13 


0.27E-03 


0.80E-08 


0.14E+00 


0.55E-04 


CPU time [s] 


- 


2.1 


36.2 


3.7 


18.0 


2.5 


12.1 


L = 18, m = 40 
















Error (N=5) 


- 


* 


0.97E-06 


0.90E-04 


0.12E-07 


0.21E-01 


0.94E-0.5 


CPU time [s] 


- 


82 


316 


74.9 


312.5 


52.3 


239.6 


L = 22, m = 8 
















Error (N=5) 


- 


* 


0.40E-06 


0.35E-04 


0.21E-07 


0.57E-02 


0.39E-05 


CPU time [s] 


- 


826 


1817 


402 


1510 


412 


1549 





networks of more complicated quantum algorithms such as Grover's database search algorithm and Shor's factoring 
algorithm These more complicated quantum algorithms are discussed in Subsections IIVI B - IIVI H. In this Section we 
discuss only the salient features of the quantum algorithms. Much more detailed information can be found in, for 
example, references S IH3, llllil HI ■ 



TABLE III: Overall comparison of different integration methods to solve the time-dependent Schrodinger equation. D = 2 
is the dimension of the Hilbert space, T is the time interval of integration, r is the time step, and A'^ denotes the number of 
iterations in the short-iterative Lanczos algorithm. An entry MP indicates that the result is numerically exact to machine 
precision. 
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A. Elementary Gates 



Hadamard Gate 



The Hadamard gate is a single-qubit gate and is often used to prepare the state of uniform superposition y\ . The 
Hadamard operation on qubit j is defined by 



'*'' = ^ 



1 1 
1 -1 



(75) 



For example 



T^,|0) = -^(|0) + |1)). 



(76) 



In terms of the elementary rotations X and Y ^ the Hadamard operation on qubit j reads 



W, = -iX^Y, 



-.y,x/ 



^X]Y, 



^Y.X]. 



(77) 



The Hadamard operation can be generalized to an arbitrary number of qubits [0| . The generalized operation is known 
as the Hadamard transform or as the Walsh-Hadamard transform. This operation consists of L Hadamard gates acting 
in parallel on L qubits. For example, for i = 2, W^2M/i|00) = (|00) + |01) + |10) + |ll))/2. The Walsh-Hadamard 
transform produces a uniform superposition of all basis states. A symbolic representation of the Walsh-Hadamard 
gate is given in Fig. QJi. 

2. Swap Gate 

The swap gate interchanges two qubits and is useful in, for example, quantum algorithms that perform Fourier 
transformations ^. In matrix notation, the SWAP operation reads 
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(78) 



showing that the SWAP operation can be decomposed into three CNOT operations. A graphical representation of 
the swap gate is shown in Fig.^). 



3. Toffoli Gate 

The Toffoli gate is a generalization of the CNOT gate in the sense that it has two control qubits and one target 
qubit |3,II3l- The target qubit flips if and only if the two control qubits Ql — Q^ — 1 (see Table llV|l . Symbolically 
the Toffoli gate is represented by a vertical line connecting two dots (control bits) and one cross (target bit), as shown 
in Fig.^. There are many ways to decompose the Toffoli gate in one- and two-qubit operations ^]. We discuss two 
examples. 

A quantum network for the first implementation is given in Fig. 01 13, 113 ■ It consists of two CNOT gates and three 
controlled phase shifts. 




FIG. 4: Quantum network for the Toffoli gate using CNOT gates and controlled phase shifts, which are part of the operations 
V and V. 




FIG. 5: Quantum network for the Toffoli gate using CNOT gates and single-qubit operation A. A {A) is a rotation by 7r/4 
(— 7r/4) about the y-axis. 

The internal operation of this quantum network is shown in Table IIVI The first two columns give all possible 
combinations of Qf and Q^. The next five columns schematically show the corresponding operations as we move 
through the network (see Fig.QJ from left to right. The columns labeled "CNOT" show only the value of the second 
qubit because the first qubit acts as a control bit and does not change. The last column summarizes the operation 
on the third qubit. From Table Hvl it immediately follows that operation V has to be constructed such that V^ flips 
Qa and that VV is equal to the identity matrix. These conditions are fulfilled by taking 
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TABLE IV: Internal operation of the quantum network, Fig. 2] representing the Toffoh gate. 



Ql 


Ql 


V 


CNOT 


V 


CNOT 


V 


Resuh 








1 


V 




V 






Ql 
VVQl = Ql 


1 







1 


V 





V 


VVQl = Ql 


1 


1 


V 







1 


V 


V^Ql = 1-Ql 



TABLE V: Internal operation of the quantum network, Fig. |^ y4 is a single-qubit rotation about the j/-axis by n/4. The gate 
N flips the qubit. 
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Taking into account that V should change only the target qubit (2) if the control qubit (1) is one, we find 
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The matrix V12 can be brought into a similar form as the controlled phase shift: 
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(81) 
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The controlled phase shift R2i{tt/2) can be implemented on the Ising model quantum computer (|20|1 by putting 
h = -J/2 and tJ = -7r/2. 

It is easy to see that the quantum network shown in Fig. ^corresponds to the following sequence of operations: 



TOFFOLI = Y3R31Y3Y2I21Y2Y3R32Y3Y2I21Y2Y3R32Y3, (82) 

where I21 = /2i(7r), Rji = Rji(Tr/2), and Rji = i?ji(— 7r/2). The sequence ()82|l can be shortened by observing that 
Y2I21Y2Y3 = Y3Y2I21Y2. This leads to 



TOFFOLI = Y3R31Y2I21Y2R32Y2I21Y2R32Y3. 



(83) 



In Fig. 13 we show a quantum network that performs the same operation as the Toffoli gate up to some known 
phase factors |3,|^. The internal operation of this network is summarized in Table Ivl The first two columns give all 
possible combinations of the control qubits Qf and Q2 ■ The next seven columns schematically show the individual 
operations as we move through the network from left to right (see Fig. (SJ. Entries in the colums labeled "CNOT" 
having a value "N" indicate that the target qubit (qubit 3) flips. The last column shows the full operation on the 
third qubit. From Tabled it follows that the first and second row of the last column are trivially satisfied. The third 
row implies that A^ = 1, which is impossible. Nevertheless, if we choose A — e^^^"/'^ we find that the condition in 
the fourth row is satisfied and that A AN A A = 2YS^Y = —2S^. This means that the target qubit will acquire an 
extra phase factor (—1) if and only if the first control qubit is one and the second control qubit is zero. Thus, the 
quantum network of Fig. [3 does not implement a Toffoli gate but performs an operation that is very similar, up to a 
phase factor that depends on the state of the control qubits [3 • 



B. Quantum Fourier Transform 



The quantum Fourier transform (QFT) is an essential subroutine in, for example, Shor's factoring algorithm, the 
order finding algorithm and phase estimation in general [3, 1^ ■ The QFT is not a new kind of Fourier transform: 
it is a particular application of the standard discrete Fourier transform. The latter is defined by x = Uy where 



{xo,...,XN-i), y = (yo,- 



,2/Af- 



_i) and Uj^k = e^'^^^^l^ I^TN . Clearly [/ is a unitary matrix. In the context of 



quantum computation, the vectors x and y are just two different representations of the basis vectors that span the 
Hilbert space of dimension iV. The QFT is defined as 



Af-l 



JV-1 
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(84) 



where the amplitudes hk are the discrete Fourier transform of the amplitudes Oj and N = 2^ where L is the number 
of qubits. Readers more familiar with traditional quantum mechanics recognize the similarity with the coordinate 
and momentum representation. Exploiting the massiveparallelism of the ideal quantum computer, the QFT can be 
carried out in 0(log2 A^) = C>{L) quantum operations |9ll3g|. 

How a quantum network can be derived for the QFT is explained in, for example, Ref. j^ilSg. In Fig.|B| we show a 
quantum network that performs a four-qubit QFT ^f7|. The blocks labeled W perform a Walsh-Hadamard transform, 
and the other blocks perform a controlled phase shift by the angle indicated. Not shown in the network is the series 
of SWAP-gates that interchange the output qubits (1,4) and (2,3) @,|3^, hence the different labehng of input and 
output lines in Fig. For the applications that we discuss later, these interchanges merely add to the operation 
count and can therefore be omitted. To sec how the network carries out the QFT (and also for the derivation of the 
network), the states \xj) and \yj) in Eq. H84|) have to be written in binary representation, for example, \xo) = |0000), 
|a;i) = |0001), and so on. 
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FIG. 6: Quantum network that performs a four-qubit quantum Fourier transform. W denotes the Walsh-Hadamard transform. 
The operations "2", "4", and "8" perform controlled phase shifts with angles 7r/2, 7r/4, and 7r/8, respectively. 



C. Finding the Period of a Periodic Function 

Assume that we are given the function j{n) = f{n + M) for n = 0, . . . , iV — 1. On a quantum computer, we can 
determine the period M as follows. We use one register of qubits to store n and another one to store f(n). As a 
first step, the quantum computer is put in the state of uniform superposition over all n. The state of the quantum 
computer can be written as 
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(85) 



ri=0 



where, in the last step, we used the periodicity of f{n). Using the Fourier representation of \n) we obtain 
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A:=0 n=0 



where L = IN/M\ denotes the largest integer L such that ML < N. In simple terms, L is the number of times the 
period AI fits into the interval [0, iV — 1]. The probability pq{M) to observe the quantum computer in the state \q) is 
given by the expectation value of the (projection) operator Q = |g)(9|. With the restriction on f{n) that /(n) = f{n') 
implies n = n' , we find 
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The results for Pq{M) in the case TV = 8 (3 qubits) are given in Table IVll From Table IVll it follows directly that the 
expectation values of the qubits are (Qf = Q2 = Q3 = 0) if the period M = 1, (Qf = Q2 — 0, Q3 — 0.5) if the period 
M = 2, (Qf = 0.5, gi = 0.375, Q| = 0.34375) if the period M = 3, and (Qf = 0, Q| = Q| = 0.5) if the period 
M = 4. Thus, in this simple case the periodicity of f(n) can be unambiguously determined from the expectation 
values of the individual qubits. 

TABLE VI: Probability pq{M) to observe the state \q) after performing the quantum Fourier transform on the perodic function 
/(n) = f{n + M) forn = 0, . . . , 7. 
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P,{M = 4) 



0.5 
0.0 
0.0 
0.0 
0.5 
0.0 
0.0 
0.0 



0.34375 
0.01451 
0.06250 
0.23549 
0.31250 
0.23549 
0.06250 
0.01451 



0.25 
0.00 
0.25 
0.00 
0.25 
0.00 
0.25 
0.00 



D. Grover's Database Search Algorithm 

Next we consider Grover's database search algorithm to find the needle in a haystack |a, Q- On a conventional 
computer, finding an item out of N elements requir es 0{ N) queries (^. Grover has shown that a quantum computer 
can find the item using only 0{\/N) attempts [M Ifl. Il0(l | . Assuming a uniform probability distribution for the needle, 
for iV = 4 the average number of queries required by a conventional algorithm is 9/4 |la.l5(]|. With Grover's quantum 
algorithm, the correct answer for this database with four items can be found in a single query |l4l llq. Grover's 
algorithm for the four-item database can be implemented on a two-qubit quantum computer. 

The key ingredient of Grover's algorithm is an operation that replaces each amplitude of the basis states in the 
superposition by two times the average amplitude minus the amplitude itself. This operation is called "inversion 
about the mean" and amplifies the amplitude of the basis state that represents the searched- for item [^ |^{ . To see 
how this works, it is useful to consider an example. Consider a database containing four items and functions gj{x), 
j — 0, ... ,3 that upon query of the database returns minus one ii x = j and plus one ii x ^ j. Let us assume that 



the item to search for corresponds to, for example, 2 (52(0) = 52(1) — 52(3) = 1 and .92(2) = — 1). Using the binary 
representation of integers, the quantum computer is in the state (up to an irrelevant phase factor as usual) 



1$) = _(|00) + 101) - |10) + 111)). 

The operator B that inverts states like (|88|l about their means reads 

/ ao 



(88) 





(89) 



Applying B to |<f>) results in B\^) — |10), that is, the correct answer. In general, for more than two qubits, more than 
one application of B is required to get the correct answer [^.Q. In this sense the two-qubit case is somewhat special. 
Implementation of this example on a two-qubit quantum computer requires a representation in terms of elementary 
rotations of the preparation and query steps and of the operation of inversion about the mean. Initially, we bring 
the quantum computer in the state |00) and then transform |00) to the state H88(l by a two-step process. First we 
use the Walsh-Hadamard transform to bring the quantum computer in the uniform superposition state: M^2W^i|00) — 
(|00) + |01) -I- |10) -I- |ll))/2, where Wj is given by Eq. H75f) . Next we apply a transformation F2 that corresponds to 
the application of 52 {x) to the uniform superposition state 
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This transformation can be implemented by first letting the system evolve in time 



(90) 



GW^2W^i|00) = 



-irrSf So 



1 



(|00) + |01) + |10) + 111)) = Tr(e~"/^|00) + e+*"/^|01) + e+'"/*|10) + e-'^/^\n)), (91) 



where the two-qubit operation G is defined by 
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(92) 



Then we apply a sequence of single-spin rotations to change the four phase factors such that we get the desired state. 
The two sequences YjXjYj and YjXjYj operating on qubit j are particularly useful for this purpose since 



r,X,y,|0)=e+-/*|0), Y,X,Y,\1) 



-iTr/4 



1), YX,Y, 



-i-n/4 
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(93) 



We find 



Y1X1Y1Y2X2Y2 
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(94) 



Thus we can construct the sequence Fj that transforms the uniform superposition state to the state that corresponds 
to gj{x): 



Fq — YiX lY 1Y2X 2Y 2G , Fi — YiX lY 1Y2X2Y 2G , F2 — YiXiY 1Y2X2Y 2G , F3 — YiXiY 1Y2X2Y 2G . 



(95) 



Finally, we need to express the operation of inversion about the mean - that is, the matrix B [see (|89|l ] - by a 
sequence of elementary operations. It is not difficult to see that B can be written as 



B = W1W2 
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0-100 
0-10 
0-1 



WIW2 = W1W2PW1W2. 



(96) 



The same approach that was used to implement 52(2;) also works for P (= —Fq) and yields 

P = -Y1X1Y1Y2X2Y2G. (97) 

The complete sequence Uj operating on |00) reads 

Uj = W1W2PW1W2FJW2W1. (98) 

Each sequence Uj can be shortened by observing that in some cases a rotation is followed by its inverse. Making use 
of the alternative representations of the Walsh-Hadamard transform Wi [see Eq. H77|l ]. the sequence for, for instance, 
j — I can be written as 

W1W2F1 = X1X1Y1X2X2Y2Y1X1Y1Y2X2Y2G = X1Y1X2Y2G. (99) 

The optimized sequences Uj read 

Uq = —X1Y1X2Y2GX1Y1X2Y2GX2X2Y2X1X1Y1, 

Ui = —XiY 1X2Y2GX1Y1X2Y 2GX2X2Y 2X1X1Y I, 
U2 = —XiYiX2Y2GXiYiX2Y2GX2X2Y2XiXiYi^ 
U3 - -X1Y1X2Y2GX1Y1X2Y2GX2X2Y2X1X1Y1. (100) 

Note that the quantum algorithms Ij 100(1 are by no means unique: various alternative expressions can be written down 
by using the algebraic properties of the Xs and Ys. Straightforward calculations show that Uj\Q) = \j). Hence the 
application of a sequence IjlOOfl to the initial state |0) brings the quantum computer in the state that corresponds to 
the searched-for item. 

E. Finding tiie Order of a Permutation 

We consider the problem of finding the order of a permutation. An experimental realization of this quantum 
algorithm on a five-qubit NMR quantum computer for the case of a permutation of four items is described in Ref . |24| | . 
The problem is defined as follows: Given a permutation P of the integers {0, 1, ..., N— 1} and an integer < y < N—1, 
the order r{y) is the smallest integer for which P^{y)y ~ y. Thus, the purpose of the quantum algorithm is to determine 
r{y), for a given y and permutation P. The theory in this section closely follows Ref. |24|. We therefore also consider 
the case TV = 4, as an example. 

The basic idea of the quantum algorithm to find the order of a permutation is to exploit the quantum parallelism 
to compute P{y), P^{y), P^{y), and P'^{y) simultaneously and filter out the power r{y) that yields P^{y)y = y. 
Denoting f{n) — P", finding r{y) is the same as finding the period of /(n), a problem that is solved by invoking the 
QFT (see Section rrvTl . 

First we consider the problem of generating a permutation of A^ = 4 items. We need two qubits to specify one of 
the four items. Using the binary representation of the integers {0, 1, 2, 3}, it is easy to see that the CNOT operation 
C21 (the right subscript denoting the control bit) corresponds to the permutation (in cycle notation) P = (0)(2)(13), 
that interchanges items 1 and 3. Likewise, C12 generates the permutation P — (0)(1)(23) and C12C21C12 is equivalent 
to the permutation P = (0)(3)(12). The remaining interchanges of two items can be generated by a combination of 
CNOT gates and NOT operations. Denoting the NOT operation on the jth qubit by Nj, we find that C21N2 = N2C21 
yields P = (02)(1)(3), G12N1 = N1G12 yields P = (2)(3)(01), and G12N2C21C12 yields P = (1)(2)(03). Using these 
elementary interchanges, we can construct any permutation. 

The quantum network that carries out the quantum algorithm to find the order of a permutation P of four items is 
shown in Fig. [7| There is a one-to-one mapping from this network onto the quantum algorithm to find the period of 
the function /(n) (see Section liVCjI . The first three qubits hold the integer n; qubits 4 and 5 hold y = 2yi +yo. The 
three Walsh-Hadamard operations change the initial state \000)\yiyo) into \uuu)\yiyo) where we use the label "m" to 
denote the uniform superposition. Then, we apply the permutations P^y, P^y, ■ ■ ■ , P'^y- In the actual implementation, 
the sequence of CNOT and NOT operations that implement the permutation P need to be replaced by Toffoli and 
CNOT gates, respectively, because (the power of) P is applied to the fourth and fifth qubit (representing the integer 
y), conditional on the state of the first three qubits. Finally, we perform a QFT on the first three qubits and consider 
the value of the first three qubits Qij Q2' ^^"^ Ql- ^^ explained before, in this simple case we can extract the period 
of the function, and hence the order of P acting on y, from the values of Qf , Q|, and Q3. 
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FIG. 7: Quantum network of a quantum algorithm to find tfie order of a permutation of four items. W and P* denote 
the Walsh-Hadamard transform and k applications of the permutation P, respectively. The operations "2" and "4" perform 
controlled phase shifts with angles 7r/2 and 7r/4, respectively. 



F. Number Factoring: Shor's Algorithm 



As another application of the QFT, we consider the problem of factoring integers, that is, Shor's algorithm. For 
the case N = 15, an experimental realization of this quantum algorithm on a seven-qubit NMR quantum computer is 
given in Ref. [23|. The theory in this section closely follows Ref. pU^ . 

The theory behind Shor's algorithm has been discussed at great length elsewhere jaSISal- Therefore, we recall 
only the basic elements of Shor's algorithm. Shor's algorithm is based on the fact that the factors p and q of an integer 
N = pq can be deduced from the period M of the function f{j) = a^modN ior j — 0, . . . , 2" — 1 where N < 2". Here 
a < A'^ is a random number that has no common factors with N. Once M has been determined, at least one factor 
of N can be found by computing the greatest common divisor (g.c.d.) of N and a^^/^ ± 1. 

Compared with the example of the previous section, the new aspect is the modular exponentiation a-'modA^. 
For A^ = 15 this calculation is almost trivial. Using the binary represention of j, we can write a^modA^ = 
a^ ■'"-1 . . . a^^^a-'^modA^ = (a^ ^"-^modA^) . . . (a^-'imodA)(a-'fmodA^)modA^, showing that we only need to im- 
plement (a2'°-"=modAf). For A^ = 15 the allowed values for a are a = 2,4,7,8,11,13,14. If we pick a = 2,7,8,13 
then a^ modA^ = 1 for all fc > 1 while for the remaining cases we have a^ modA^ = 1 for all fc > 0. Hence for all 
a we need to "calculate" ImodA and amodA^ and for a = 2, 7, 8, 13 we in addition have to calculate a^modA^ and 
a^modA". Thus, for A^ — 15, only two (not four) qubits are sufficient to obtain the period of f{j) = a^iaodN p3|. As 
a matter of fact, this analysis provides enough information to deduce the factors of A^ = 15 using Shor's procedure 
so that no further computation is necessary. Nontrivial quantum operations are required if we decide to use three (or 
more) qubits to determine the period of /(j) — a^uiodN |20|- Following Ref. |23> we consider a seven-qubit quantum 
computer with four qubits to hold f(j) and three qubits to perform the QFT. 

In essence, the quantum network for the Shor algorithm is the same as the one shown in Fig. {7\ (and therefore not 
shown) with the permutations (two qubits) replaced by modular exponentiation (four qubits). The quantum networks 
to compute a-'modlS for j = 0, . . . , 7 and a fixed input a are easy to construct. For example, consider the case a = 11 = 
|1011). If j is odd, then ll-'modl5 ~ 11 and the network should leave |1011) unchanged. Otherwise, ll-'modl5 = 1, 
and hence it should return |0001) (in this case M = 2 and g.c.d.(A^, o^/^i 1) = {g.c.d.(15, 10), g.c.d.(15, 12)} = {5, 3}, 
showing that there is no need to perform a quantum computation). The network for this operation consists of two 
CNOT gates that have as control qubit, the same least-significant qubit of the three qubits that are input to QFT. 
The sequence of CNOT and Toffoli gates that performs similar operations for the other cases can be found in the 
same manner. 
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FIG. 8: Quantum network of a three-input adder, as described in Ref. |101|| . The algorithm performs a quantum Fourier 
transform of register 3 (qubits 9 to 12), adds the content of register 2 (qubits 5 to 8) and the content of register 1 (qubits 1 
to 4), followed by a QFT on register 3 to yield the final answer (register 1 + register 2 + register 3 modl6). W denotes the 
Walsh-Hadamard transform. The operations "±1," "±2," "±4," and "±8" perform controlled phase shifts with angles ±7r, 
±7r/2, ±7r/4 and ±7r/8, respectively. Squares that touch each other indicate operations that can be done simultaneously. 



G. A Three-Input Adder 



This subsection gives another illustration of the use of the QFT: a quantum al gorit hm to add the content of three 
four-qubit registers. This example is taken from the Ph.D. thesis of S. Bettelli |lOl| . The quantum network of the 
complete circuit is shown in Fig.|Sl The modular structure of this approach is clear. Note that with respect to the QFT 
network of Fig. |21 both the labeling of qubits and the order of operations have been reversed. The former is merely 
a change of notation and the latter allowed because quantum algorithms are reversible (unitary transformations) by 
construction Q. The basis idea of this algorithm is to use the QFT to first transfer the information in a register 
to the phase factors of the amplitudes and then use controlled phase shifts to add information from the two other 
registers and finally QFT back to the original representation. Note that this quantum network differs considerably 
from the one described in Ref. [^ and is also easier to implement. 



H. Number Partitioning 



As a final exa mple , we discuss a quantum algorithm to count the number of solutions of the number partitioning 
problem (NPP) |102| . The NPP is defined as follows: Does there exist a partitioning of the set A — {ai, ..., a„} of n 
positive integers gj into two disjoint sets Ai and A2 — A — Ai such that the difference of the sum of the elements of 
Ai and the sum of the elements of A2 is either zero (if the sum of all elements of A is even) or one (if the sum of all 
elements of A is odd)? The following simple example may be useful to understand the problem. li A = {1,2,3,4}, 
the answer to the NPP is yes because for Ai = {1, 4} and A2 = {2, 3}, the sum of the elements of Ai and A2 are both 
equal to five. In this case there are two solutions because we can interchange Ai and A2. If A = {1, 1, 1, 4} the answer 
is again yes, as there is one solution, namely Ai = {1,1,1} and A2 = {4}. The difference of the sum of the elements 
of Ai and A2 is equal to one, and that is all right because the sum of all elements of A is odd. li A = {2, 2, 2, 4}, there 
is no solution to the NPP. 

The quantum network for the NPP quantum algorithm for the case that the sum of four integers to be partioned is 
maximally equal to 16 is shown in Fig.lHl The basic idea behind this quantum algorithm is a number of transformations 
that reduce the c ounti ng of the number of solutions of the NPP to finding the number of zero eigenvalues of a spin- 
1/2 Hamiltonian |l02| . The largest part (in terms of the number of operations) of the quantum algorithm is a 
combination of gates that transforms the state of the quantum computer such that the number of solutions n^ i s a 
physically measureable quantity (the expectation values of the 15th spin). For n = 4, we have ris = 16^/Qf^ |l02j . 
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FIG. 9: Block diagram of the quantum algorithm that solves the number partitioning problem. The first four qubits are used 
to represent the sum of the four integers to be partitioned. Qubits 5 to 8 are used to determine the number of solutions of 
the number partitioning problem. The remaining seven qubits are used to relate the number of solutions Us to a physically 
measurable quantity: the expectation value of the 15th qubit. The unitary transformation U prepares the uniform superposition 
of the first eight qubits, U is the inverse of U , and the combination of INVERT and AND gates sets the 15th qubit to one if 
and only if the first eight qubits are all one. 

V. SIMULATION OF IDEAL QUANTUM COMPUTERS 

On the most abstract level, a quantum algorithm that runs on an ideal quantum computer performs one unitary 
transformation {D x D matrix) of the state vector (of length D = 2^ , where L is the number of qubits). On a 
conventional computer, this calculation would take 0{D^) arithmetic operations. As in the case of programming a 
conventional computer, it is extremely difficult to write down this one-step operation explicitly. Usually an algorithm 
consists of many steps, each step being a simple construct that easily translates into the elementary operations of the 
processing units. 

As explained in Section III CI for a quantum computer the set of single-qubit rotations and the CNOT gate or any 
other set of gates that can be used for universal quantum computation will qualify as elementary operations. Thus, 
a first task is to decompose a given quantum algorithm (unitary matrix) into a sequence of unitary matrices that act 
on one or two qubits. Recall that each of these matrices is very sparse (see Section Fill All . Multiplying a very sparse 
matrix and a vector can be done very efficiently on a conventional computer, typically in 0{D) arithmetic operations. 
Hence the remaining important question is whether a quantum algorithm can be broken down such that the number 
of elementary operations is much less than D. The general answer to this question is no |2|. However, for all quantum 
algorithms reviewed in this chapter, the answer is affirmative. Finding the shortest decomposition of a given unitary 
matrix in terms of a fixed set of sparse unitary matrices is a difficult optimization problem. 

Assuming that we have a representation of a quantum algorithm in terms of the elementary set of sparse unitary 
matrices, simulating this algorithm on a conventional computer is conceptually simple. The basic techniques are 
reviewed in Section fHI Al Software that performs these calculations is called a quantum computer gate-level simulator. 
It is an essential tool to validate the gate-level design of a quantum algorithm on an ideal quantum computer. Examples 
of such simulators are discussed in Section IVlII 



VI. SIMULATION OF PHYSICAL MODELS OF QUANTUM COMPUTERS 

The qubits on an ideal quantum computer are ideal two-state quantum systems. Therefore, the operation of an ideal 
quantum computer does not depend on the intrinsic dynamics of its qubits. However, a physically realizable quantum 
computer is a many-body system of qubits. The quantum dynamics of these qubits are at the heart of real quantum 
computation. Conceptually, the main differences between simulating quantum computers on a quantum-gate level 
(see Section^ and physical models of quantum computers are that the dynamics of the latter are specified in terms of 
a time-dependent Hamiltonian (not in terms of unitary matrices) and that the physical quantum computer is a single 
physical system (not a collection of abstract mathematical entities that can be controlled with infinite precision) . 



T here are many prop osals for building quantum computer hardware [T^, Il03l Il04 llOSl Il06l llOTl llOSl Il09l IllOl 

lllll III2I lll.li Ill4 'll5l , and the criteria that a physical system should satisfy to qualify for quantum computation 
have been examined 116 , 117 , 118, 119] . Candidate technologies for building quantum gates include ion traps, ca vity 
QED, Josephson junctions, and nuclear mjimetic resonanc e (NMR) technolog y [1 11^ 11^ 111 111 11^ Im [^ 

pl ^El E l p M E llll lSliaiaiMlMiaiMIMGMlM 

I121I I122. Il2,'l [124I Il2f)l Il2fil |. With the exception of NMR techniques, none of these technologies has demonstrated 
nontrivial (that is, more than one CNOT operation on an arbitrary linear superposition) quantum computation. In 
this respect, the field of simulating quantum computers is mu ch more developed: it is relatively easyto simulate 
realistic physical models of NMR quantum computers |93lll27j | or a pair of Josephson junction qubits |l28lll29| |. 

In some cases (for instance, electron spins, nuclear spins), there is a one-to-one mapping from the mathematical 
object of a qubit to a dynamical variable in the quantum mechanical description. In other cases (for example, 
Josephson junct ions), thi s mapping is not self-evident and requires a very careful analysis of the dynamics of the 
physical system |l28l Il2 9j . Mathematically, nothing prevents us from using three- or even many- level systems as basic 
units for quantum computation. However, inventing useful quantum algorithms that use two-state systems as qubits 
already seems to be a major tour de force. Therefore, it is not clear whether working with three- or more-valued 
qubits instead of two-state systems will bring more than additional complications. Thus, in the following we take the 
point of view that the quantum computation will be carried out with physical systems that are described by (but not 
necessarily are) two-state quantum systems. 

The physical system defined by Eq. H21(l is sufficiently general to serve as a physical model for a generic quantum 
computer at zero temperature. For instance, Eq. H21|l includes the simplest (Ising) model of a universal quantum 
computer |5lll53| 

L 3 

i,j=l j = l a=x,y,z 

More specific candidate hardware realizations of Eg. (12 111 include linear arrays of quantum dots |lllj| , Josephson 
junctions ^3: ^'^^ NMR systems [Ij, [I^ [la, llTl [251 [2^ [271 [23 |. An approximate model for the linear arrays of 
quantum dots reads 

L L L 

H(t) = - y^ i?,5|5|+i - ^ h^^{t)S^ +EoyZ P,{t)Sl (102) 



i=i j=i i=i 



where Ej = Eq {Ej = 2Eq) when j is odd (even) and /ij(t) and Pj{t) are external control pa rameters Jl 1 ij ■ 
Projection of the Josephson-junction model onto a subspace of two states per qubit yields |l04lll2Cl| 



L L L 

Hit) = -2Ej{t) J2 Sp]+^ -EjY^S^-Y, h^it)S], (103) 

where the energy of the Josephson tunneling is represented by Ej and Ej (t) denotes the energy associated with the 
inductive coupling between the qubits |104 , .1201 1 . Here ft.| (t) and Ej (t) may be controlled externally. 

For nuclei with spin quantum number S = 1/2, the Hamiltonian that describes the int eracting s pin system is of 
the form l|21|l 130] . NMR uses radio- frequency electromagnetic pulses to rotate the spins |l30l Il3ll| . By tuning the 
radio frequency of the pulse to the precession (Larmor) frequency of a particular spin, the power of the applied pulse 
(= intensity times duration) controls how much the spin will rotate. The axis of the rotation is determined by the 
direction of the applied RF-field. By selecting the appropriate radio-frequency pulses, arbitrary single-spin rotations 
can be carried out. In other words, using radio-frequency pulses we can perform any single-qubit operation. The 
simplest model for the interaction of the spins with the external magnetic fields reads 

h'^it) = ^ + h"; sin(27r/j"t + (^p, (104) 

where /i" and /i" represent the static magnetic field and radio-frequency field acting on the jth spin, respectively. 
The frequency and phase of the periodic field are denoted by /" and </?" . The static field ft." fixes the computational 
basis. It is convention to choose M = hH — Q and M 7^ so that the direction of the static field corresponds to the 



z-axis of the spins. Under fairly general conditions and to a very good approximation, the nuclear spin Haniiltonian 
(containing exchange interactions, dipole-dipole interactions, and so on) reduces to the universal quantum computer 
model H20|l in which there are interactions only between the z-components of the spin operators. All NMR quantum 
computer experiments to date have been interpreted with this model. Therefore, it makes sense to use model (|2U|I as 
a starting point for simulating physical realizations of quantum computers. 

A. NMR-like Quantum Computer 

In this section, we illustrate the difference between simulating an ideal, c ompu ter science-type quantum computer 
and a more realistic, physical model of quantum computer hardware |93l Il27j . We limit our presentation to the 
implementation of the CNOT gate and Grover's algorithm on a two-qubit, NMR-like quantum computer. A more 
extensive discussion, as well as many other examples, can be found in Ref. |l27j . All simulations have been carried 
out using the quantum computer emulator software |132i |. The examples given in this section are included in the 
software distribution of the quantum computer emulator software |l32l | . 

As a prototype quantum computer model, we will take the Haniiltonian for the two nuclear spins of the ^H and "'^^C 
atoms in a carbon-13 labeled chloroform molecule that has been used in NMR-quantum computer experiments |lfiLll7l |. 
The strong external magnetic field in the z-direction defines the computational basis. If the spin is aligned along the 
direction of the field, the state of the qubit is |0); if the spin points in the other direction, the state of the qubit is 
|1). In the absence of interactions with other degrees of freedom, the Hamiltonian of this spin-1/2 system reads 

Hnmr = —Ji,2^i^2 ^ ^i'5'i ~ h2S2, (105) 

where h\/2n k, SOOMHz, ft,|/27r « 125MHz, and J = Jf 2/27r w — 215IIz 16]. In our numerical work, we use the 
rescaled model parameters 

J = -0.43xl0"^ hl^l, hl^l/A. (106) 

The ratio 7 — /i2/^i — 1/4 expresses the difference in the gyromagnetic ratio of the nuclear spin of the ^H and ^^C 
atom. 

1. Single-Qubit Operations 

In NMR experiments, it is impossible to shield a particular spin from the sinusoidal field. An application of a 
sinusoidal field not only affects the state of the resonant spin but also changes the state of the other spins (unless they 
are both perfectly aligned along the z-axis). An analytical, quantitative analysis of this simple-looking many-body 
problem is rather difficult. The values of the model parameters Hl()6(l suggest that the interaction between the spins 
will have a negligible impact on the time evolution of the spins during application of the sinusoidal pulse if the duration 
of the pulse is much shorter than 1/ J (we use units such that tJ is dimensionless). Thus, as far as the single-qubit 
operations are concerned, we may neglect the interaction betwee n the two spins (which is also confirmed by numerical 
simulation of l|105|l . see later). In this section, we closely follow jl27| . 

We consider the two-spin system described by the time-dependent Schrodinger equation 

i|-|$(t)) = - hlS^ + ^S^ + hl{S^sinujt + SfcosLut) + h'^{S^smujt + S^cosLut) |$(i)), (107) 

for tw o int eracting spins in a static and a rotating sinusoidal field. As usual, it is convenient to work in a rotating 
frame [ipl- Substituting |$(t)) = e'^'^^^'^+'^^^\-^{t)), we obtain 

'§l\^(i)) = ~[iK~^)st + m-^)si + Ksy + ^sy] |*(t)). (los) 

Our aim is to determine the conditions under which we can rotate spin 1 by an angle (pi without affecting the state 
of spin 2. First we choose 

u; = hl, (109) 



that is, the frequency of the sinusoidal field is tuned to the resonance frequency of spin 1. Then (|108|l can easily be 
integrated. The result is 

|$(t)) ^e«t'»i(Si+sj)e^*^Tsre**S2-^i2|$(0)), (110) 

where v„m = {0,h%^,h^ — h^). The third factor in (|ll()|l rotates spin 2 about the vector Vi2. This factor can be 
expressed as 



e-— = I 1 "■ ^= + 'I-.!-' -r,, " ,= :"A '^"^T^. 0") 



m n in 



and we see that the sinusoidal field will not change the state of spin 2 if and only if the duration ti of the pulse 
satisfies 



ii|vi2| = ti V(M - hi? + m? = 4^"i, (112) 

where ni is a positive integer. The second factor in H11Q(I is a special case of Hlll|) . Putting 

tihf = y'l, (113) 

the second factor in (|110|l will rotate spin 1 by ipi about the y-axis. Therefore, if conditions H109|l . (|112|l . and H113() 
are satisfied, we can rotate spin 1 by ipi without affecting the state of spin 2, independent of the physical realization 
of the quan tum computer. Combining these conditions yields the so-called 27rfc method for suppressing nonresonant 
effects nMllsl: 



IM(1~7)I 



^16nl - (^i/7r)2. (114) 



If we would use the conditions Hl()9(l . 1113|l and l|114|l to determine the parameters ui, hf and ti of the radio- 
frequency pulse, the first factor in Hll()|l can still generate a phase shift. This phase shift clearly depends on the 
state of the spins. Although it drops out in the expression of the expectation value of single qubits, for quantum 
computation purposes it has to be taken into account (as is confirmed by numerical calculations ,127. |). Adding the 
condition 

tihl^Airki, (115) 

where fci is a positive integer (/if > by definition), the first factor in 1110() is always equal to one. A last constraint 
on the choice of the pulse parameters comes from the fact that 

h^^jK , h^^jh^ ■ a^x,y,z. (116) 

Without loss of generality, we may assume that < 7 < 1. 

Using conditions (|109() . H112() . H113() . H115() . and (|116|l and reversing the role of spin 1 and spin 2, we obtain 

where fci, ^2, ni, and ^2 are positive integers. The angles of rotation about the y-axis can be chosen such that 
< V^i < 27r and < ip2 < 27r. Of course, similar expressions hold for rotations about the x-axis. 

In general, (|117|l has no solution, but a good approximate solution may be obtained if 7 is a rational number and fci 
and ^2 are large. Let 7 = N/AI where N and M are integers satisfying < A^ < M. It follows that the representation 
fci — kMN'^ and fc2 — kNM"^ will generate sufficiently accurate solutions of (|117(l if the integer k is chosen such that 
2kNM{M — A^) ^ 1. In terms of k, N, and M, the relevant physical quantities are then given by 



We have derived conditions (|118|l under the assumption of ideal sinusoidal RF p ulses. In a n experiment, there is no 
such limitation: the sinusoidal fields may be modulated by almost any waveform |l3(l Il34| . However, the fact that 
in quantum computer applications it is necessary to use single-spin pulses that do not change the state of the other 
spins remains. For general pulses, finding the form of the pulse that rotates spin 1 such that the state of spin 2 is not 
affected is a complicated nonlinear optimization problem J1271 Il35j . 

To summarize: If conditions 1)1091) . H112|) . H113|l . and (|115|l are satisfied, we can rotate spin 1 by ipi without affecting 
the state of spin 2 and without introducing unwanted phase shifts. In numerical experiments, Ea. (|118() may be used 
to determine the duration of the sinusoidal pulses. These sinusoidal pulses will then be optimized in the sense that a 
pulse that rotates spin 1 (2) will change the state of spin 2 (1) only shghtly if k satisfies 2kNM{M — iV) ^ 1. 

2. Two-Qubit Operations 

In Section ITlFl we implemented the CNOT sequence (EHJ using the Ising model H = -JSfSi - h{Sf + 5|). The 
implementation of the CNOT operation using Eq. 1)105(1 requires additional steps to account for the fact that the two 
nuclear spins feel different static fields. The additional rotations are 

CNOT = y2e^'''^'''i^'''^ie^*'^'''^^''^^2g-«'^-f^NMHy2 _ Y2e^^'^^^^^^''^^e~'^'^^'^^~'^^^^Y2Y2e^'^'^^"^"^Y2 (119) 

where we used the fact that Y2Y2 = 1. The extra phase shifts in H119|l can be expressed in terms of single-qubit 
operations. The identities 

^~tr{hl-h)s', ^ YiX[Yi = XiY{Xi, e-^r(h'2-h)s'2 = Y2X'^2, (120) 

define the single-spin rotations X[, Y(, and Xj. 

In the case of Grover's database search algorithm, the representation of G in terms of the time evolution of H105|l 
reads 

G = e-*^'5i'S'2 ^ ^-^ThfSf^-iThisi^-trHMMR ^ Y2X2Y2YiX'^Yie^""'^'"'^ , (121) 

where r = — tt/ J. This choice of t also fixes the angles of the rotations and all parameters of the operations X'^ and 

Equation H120|l suggests that there are many different, logically equivalent sequences that implement the CNOT 
gate on an NMR-like quantum computer. We have chosen to limit ourselves to the respresentations 

CNOTi = YiX[YiX!,Y2l'Y2, CNOT2 = yiX^X^FiFa/'Fa, CNOT3 = XiF/X^FaXi/'ya, (122) 

where we introduced the symbol /' to represent the time evolution ^-^'^Hnme. ^ffilh r = —tt/J. 

On an ideal quantum computer, there is no difference between the logical and physical computer and the sequences 
(I122|l give identical results. However, on a physical quantum computer such as the NMR-like quantum computer (|l()5(l . 
this is not the case. On a physically realizable NMR-like quantum computer X1X2 ^ -^2^1 unless 2kN M [M — N) ^ 1 
and (|118|l are satisfied exactly. Next we use the sequences (I122|l to demonstrate that this unpleasant feature of physical 
quantum computers may give rise to large systematic errors. 

3. Simulation Results 

The model parameters for the rotating sinusoidal fields are determined according to the theory outlined previously. 
We use the integer k to compute all free parameters and label the results of the quantum computer calculation by the 
subscript s = 2kMN'^. For reference we present the set of parameters corresponding to s = 8 (/c = 1) in Table IVIll 
Multiplying s (the duration of the sinusoidal pulse) with the unit of time (2 ns for the case at hand) shows that in our 
simulations, single-qubit operations are implemented by using short pulses that are, in NMR terminology, nonselective 



TABLE VII: Model parameters of single-qubit operations on an NMR-like quantum computer using rotating sinusoidal fields 
for the case {k — 1, N — 1, AI — 4); see l|118|l . Parameters of model l|l()7|l that do not appear in this table are zero, except 
for the interaction J — —0.43 x 10~^ and the constant magnetic fields hi — 1 and /12 — 0.25. The time-dependent Schrodinger 
equation is solved using a time step S/2n — O.Of . 
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TABLE VIII; Expectation values of the two qubits as obtained by performing a sequence of five CNOT operations on an NMR- 
like quantum computer. The initial states |10), |01), |11), and \singlet) = (|01) — |10))/\/2 have been prepared by starting 
from the state |00) and performing exact rotations of the spins. The CNOT operations on the singlet state are followed by a 
7r/2 rotation of spin 1 to yield a nonzero value of qubit 1. The time s — r /2-k — 2kMN'^ determines the duration and strength 
of the sinusoidal pulses through relations H118|l . see Table lyTTI for the example of the case s = 8. The CNOT operation itself 
was implemented by applying the CNOT sequence given by (11221 . On an ideal quantum computer, CNOT^ is the identity 
operation. For s = 256, all results are exact within an error of 0.01. 
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and hard. In contrast to the analytical treatment given in Section I VI A II in all our simulations the interaction J is 
nonzero. 

The two-qubit operation /' can be implemented by letting the system evolve in time according to Hamiltonian 
Hnmr, given by (|105|l . /' is the same for both an ideal or NMR-like quantum computer. Note that the condition 
tJ — — TT yields T/27r = 1162790.6977, a fairly large number [compared to ft.f = 1, see (|105|) ]. Al so n ote the digits 
after the decimal point: this accuracy is necessary for correct operation of the quantum computer |l27l |. 

As a first check we execute all sequences on an implementation of the ideal quantum computer and confirm that 
they give the exact answers (results not shown) . It is also necessary to rule out that the numerical results depend on 
the time step S used to solve the time-dependent Schrodinger equation. The numerical error of the product formula 
used by QCE is proportional to 6^ 62). It goes down by a factor of about one 100 if we reduce the time step by a 
factor of 10. Within the two-digit accuracy used to present our data, there is no difference between the results for 
S = 0.01 and S = 0.001. 



TABLE IX: Expectation values of the two qubits as obtained by running Grover's database search algorithm on an NMR-like 
quantum computer. The time s — r/2Tv = 2kMN^ determines the duration and strength of the sinusoidal pulses through 
relations (I118L see Table IVm for the example of the case s = 8. Within two-digit accuracy, all results for s — 256 are exact. 
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Item position Qf Ql Ql Q\ Q\ Qi Qf Ql Ql Q\ 
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In Table IVIIII we present simulation results for CNOT acting on one of the basis states and Y\ CNOT^ acting on 
a singlet state, using the three logically equivalent but physically different implementations CNOTi, CNOT2, and 
CNOT3 [see Ea. H122|l ]. It is clear that some of the least accurate implementations (s = 8) do not reproduce the correct 
answers if the input corresponds to one of the four basis states. Moreover, if the operations act on the exact singlet 
state, the results strongly depend on the CNOT implementation if s < 32. In agreement with the theoretical analysis 
given previously, the exact results are recovered if s is sufficiently large. On the time scale set by J, the pulses used 
in these simulations are so short that the presence of a nonzero J has a negligible effect on the single-qubit pulses. 
These simulations also demonstrate that in order for a quantum algorithm to work properly, it is not sufficient to 
show that it correctly operates on the basis states. 

In contrast to computation in the classical framework, quantum computation can make use of entangled states. At 
the point where the quantum algorithm actually uses an entangled state, the quantum algorithm is most sensitive 
to (accumulated) phase errors. As another illustration of this phenomenon, we present in Tabic IIXI some typical 
results obtained by executing Grover's database search algorithm. We used the same NMR-like quantum computer 
as for the CNOT calculations. We conclude that reasonably good answers are obtained if s > 32, in concert with our 
observations for the CNOT operation. 

B. Decoherence 

An important topic that we have not discussed so far is the effect of the interaction of the quantum computer 
with its environment (dissipation, decoherence). The general belief seems to be that engineers will be able to cope 
with systematic errors (sec Section I VI A 31 for examples) due to the physics of the qubits and that problems due 
to decoherence will be the main stumblin g bl ock for buildii ig us eful quan tum co mpute rs [3|- As most theoretical 
work [iMIIMIIMIlMIIiSIIilllilllilllSIIilllM is based on approximations, 

the validity of which needs to be established, computer simulation of physical models that include decoherence may 
be of great value to gain insight into this challenging problem. Dissipation cannot be treated within the context of the 
simulation models that we have covered in this review. Instead of solving the time-dependent Schrodinger equation, 
we have to solve the equations of motion of the full density matrix of an interacting many-body system. Although 
still feasible for a small number of qubits L, the computation time now scales with 2^^ instead of with 2^. 

In the absence of dissipation, it is straightforward to incorporate into the class of simulation models physical 
processes that lead to loss of phase information during the time evolution. In Section IlIIFI we have given some hints 
as how this can be done. Needless to say, performing such simulations is costly in terms of CPU time but the pay-off 
may be significant: the simple, highly idealized uncorrelated random processes that are being used to analyze error 
correction and fault-tolerant quantum computing are very far from being physically realistic J3!tI | . Quantum error 
correction schemes that work well on an ideal quantum computer require many extra qubits and many additional 
operations to detect and correct errors |35j . The systematic errors discussed previously are not included in the current 
model J2J of quantum error correction and fault tolerant computing. Our simulation results for a most simple NMR- 
like quantum computer demonstrate that systematic errors are hard to avoid, in particular if the number of qubits 
increases (which is a basic requirement for fault-tolerant quantum computation) . It would be interesting to simulate a 
two-qubit quantum computer that interacts with other spins and analyze the effect of decoherence on, for example, the 
CNOT operation. Such simulations are definitely within reach but, to our knowledge, no results have been reported. 



VII. QUANTUM COMPUTER SIMULATORS 
A. RevieAV 

In this section, we review software tools that run on conventional computers and can be used to study various 
aspects of quantum computation. The term "quantum computer simulator" is used in a broad sense: not all the 
software discussed in this section actually simulates a quantum computer. In recent years many quantum computer 
simulators have been developed. The level of complexity of the simulations they can perform varies considerably: 
some deal with both the quantum hardware and s oftwa re while others focus on quantum computer algorithms and 
software. An early detailed survey is given in Ref. [l5J|. In the survey that follows, we confine ourselves to software 
that was accessible via the Web at the time of writing. 

Because quantum computer hardware is not readily available, quantum computer simulators are valuable tools 
to develop and test quantum algorithms and to simulate physical models for the hardware implementation of a 
quantum processor. Simulation is an essential part of the design of conventional digital circuits and microprocessors 
in particular. It is hardly conceivable that it will be possible to construct a real quantum computer without using 
simulation tools to validate its design and operation. 

There is a very important difference in simulating conventional microprocessors and quantum processors. In conven- 
tional digital circuits, the internal working of each logical circuit is irrelevant for the logical operation of the processor. 
However, in a quantum computer the internal quantum dynamics of each elementary building block is a key ingredient 
of the quantum computer itself. Therefore, in the end the physics of the elementary units that make up the quantum 
computer have to be incorporated into the simulation model of the quantum computer. It is not sufficient to model 
a quantum computer in terms of hypothetical qubits that evolve according to the rules of a mathematical model of a 
hypothetically ideal quantum computer. Including models for noise and errors also is no substitute for the physics of 
the qubits that comes into play when one wants to make contact to a real physical system. 

Quantum computer simulators come in different flavors and have different levels of complexity. A first group of 
software deals with programming languages for quantum computers. Programming languages express the semantics 
of a computation in an abstract manner and automatically generate a sequence of elementary operations to control 
the computer. An overview of quantum programming languages is given in Table H A second group of simula- 
tors comprises the quantum compilers, as summarized in Table IXII A quantum compiler takes as input a unitary 
transformation and returns a sequence of elementary one-qubit and two-quibit operations that performs the desired 
quantum operation. Quantum circuit simulators, or gate-level simulators, form the third group of simulators. On an 
abstract level, quantum computation on an ideal quantum computer amounts to performing unitary transformations 
on a complex-valued (state) vector. A conventional computer can perform these operations equally well, provided 
there is enough memory to store all the numbers of the vector. That is exactly what quantum circuit simulators do: 
they provide a software environment to simulate ideal quantum computers. A summary of quantum circuit simula- 



tors is given in Table IXIH The fourth group, collected in Table IXIIIl consists of software that uses time-dependent 
Hamiltonians to implement the unitary transformations on the qubits. These simulators make it possible to emulate 
various hardware designs of quantum computers, that is, they simulate models for physical realizations of quantum 
computers. Simulators of this group can also be used as gate-level simulators. Finally, Table IXTVl describes a number 
of purely pedagogical software products. 

B. Example: Quantum Computer Emulator 

In this section, we briefly discuss some features of QCE (see Table IXIlll . a software tool that simulates ideal 
quantum computers and also emulates physical realizations of quantum computers. The QCE is freely distributed 
as a self-installing executable, containing the program, documentation and many examples of quantum algorithms, 
including all the quantum algorithms discussed in this review. The file help. htm _155j contains information about 
how to install and how to start the QCE. The QCE runs in a Windows environment. It consists of a simulator of a 
generic, general-purpose quantum computer, a graphical user interface, and a real-time visualization module. 

QCE's simulation engine is built on the principles reviewed in this chapter. The engine simulates the physical 
processes that govern the operation of the hardware quantum processor, strictly according to the laws of quantum 
mechanics. It solves the time-dependent Schrodinger equation H23|) by a Suzuki product-formula (see Section Fill E|) 
in terms of elementary unitary operations. For all practical purposes, the results obtained by this technique are 
indistinguishable from the exact solution of the time-dependent Schrodinger equation. The graphical user interface is 
used to control the simulator, to define the hardware of the quantum computer, and to debug and execute quantum 
algorithms. Using the graphical user interface requires no skills other than the basic ones needed to run a standard 



MS-Windows applications. The current version of QCE (8.1.1) simulates quantum computers with a maximum of 16 
qubits. 

QCE can be used to validate designs of physically realizable quantum processors. QCE is also an interactive 
educational tool to learn about quantum computers and quantum algorithms. As an illustration, we give a short 
exposition of the implementation of the three-input adder (see Fig.lSJ on an ideal quantum computer and of Grover's 
database search algorithm on an NMR-likc quantum computer. Other examples are given in R.efs. |,93> .102 . ■127> .Ififil l . 

1. General Aspects 

As described in Section ll VI a quantum algorithm for a quantum computer model (j21|l consists of a sequence of ele- 
mentary operations that change the state |<i>) of the quantum processor according to the time-dependent Schrodinger 
equation (|23|l . namely, by (a product of) unitary transformations. We call these elementary operations microinstruc- 
tions in the sequel. They do not play exactly the same role as microinstructions in digital processors. They merely 
represent the smallest units of operation of the quantum processor. The action of a microinstruction on the state |$) 
of the quantum processor is defined by specifying how long it acts (that is, the time interval it is active) and the values 
of all the Js and /is appearing in H{t) H21() . the model Hamiltonian of the quantum computer. A microinstruction 
transforms the input state |$(i)) into the output state \^{t + r)), where t denotes the time interval during which the 
microinstruction is active. During this time interval the only time-dependence of H{t) is through the time-dependence 
of the external fields on the spins. Microinstructions completely specify the particular (ideal or physical) realization 
of the quantum computer that one wants to emulate. Quantum algorithms are translated into quantum programs 
that are written as a sequence of microinstructions. The graphical user interface of the QCE has been developed to 
facilitate the specification of the microinstructions (to model the quantum computer ha rdwa re) and the execution of 
quantum programs. A detailed exposition of how to use the QCE can be found in Ref. |l57j |. 

2. Three-Input Adder on an Ideal Quantum Computer 

A QCE implementation of the three-input adder described in Section llVCl is shown in Fig.^| The left panel shows 
the microinstruction set "adder" and the three panels on the right show the quantum programs "H-l-l-l," "1-1-2-1-3," 
and "9-1-9-1-9." These can be found in the quantum program directory "adder." The microinstruction set contains all 
microinstructions that are needed to execute the quantum programs "l+H-l," "l-|-2-|-3," and "9-1-9+9" on an ideal 
12-qubit quantum computer. The results (in binary notation) of the examples (l-)-l-|-l, H-2-1-3, and 9-1-9-1-9) can 
be read off from the values of qubits 9 to 12 at the bottom of the quantum programs. Qubit 12 corresponds to the 
least significant bit. The numerical values of the qubits appear when the mouse moves over the bottom region of the 
quantum program window (red or dark gray corresponds to 1, green or light gray to 0, and greenish brown or middle 
gray to 0.5). The graphics area in Fig. ^] is too small to show all microinstructions (on the computer the scroll bars 
allow the user to open/edit all microinstructions). The same holds for the three quantum programs. This example 
suggests that programming more complicated quantum algorithms like this one should not be done by hand; in fact, 
the microinstructions and quantum programs for these examples have been generated by another computer program. 
Some of the microinstructions may look rather complicated, but that is a little misleading: Whenever it is logically 
allowed to perform operations simultaneously (see Fig. Il()|l , these operations have been put into one microinstruction. 

3. Grover's Algorithm on an NMR-like Quantum Computer 

The QCE software comes with examples of quantum programs that perform Grover's database search algorithm 
with four items on ideal and NMR-like quantum computers (|105|) . In this application, the quantum computer has two 
qubits. In Fig. ^2 we show the QCE window after loading the microinstruction set "grov-rf-k032" and the quantum 
programs "grovO", . . ., "grov3" from the QCE program directory "grover2." The quantum programs "grovO", . . ., 
"grov3" implement the Grover algorithms for which the searched-for item corresponds to item 0, 1, 2 or 3, respectively. 
For printing purposes, quantum program "grov3" has been omitted from Fig. 1111 The microinstruction set "grov- 
rf-k032" corresponds to the parameter set for the rotating sinusoidal fields labeled with s = 32 in Table Hxl In the 
microinstruction set "grov-rf-k032," quantum program "h_cor_l" and micro instruction "h_nmr" are used to construct 
operation G [see (I121|l ]. Running the four quantum programs yields the values for Q\ and Qf given in colums eight 
and nine of Table Hxl In contrast to the example of the three-input adder, we use the inverse binary notation: qubit 
1 corresponds to the least significant bit. 
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FIG. 10: QCE implementation ol the quantum network of Fig.|H|for adding three quantum registers of four qubits each on an 
ideal 12-qubit quantum computer. Quantum programs (from left to right) compute 1+2+3, 1+1 + 1, and 9+9+9 mod 16. Not 
shown are the microinstructions that set the initial values of the three four-qubit registers. 



QCE has an option to visualize the time evolution of the state of the quantum computer in terms of arrows 
representing the expectation values of the qubits. In Fig. ^2 this option was used to visualize the outcome of the 
program "grov2." Microinstruction sets corresponding to the parameter sets labeled with s — 8, s — 16, s — 64, and 
s = 256 are also provided with the QCE software. Running the programs "grovO", . . ., "grov3" with these instruction 
sets demonstrates that the exact results are recovered if s is sufficiently large (see colums two and three in Table Hxll . 



VIII. SUMMARY AND OUTLOOK 



The simulation methods reviewed in this chapter have the potential to faithfully simulate quantum computers of 
20 - 30 qubits, considerably more than what is projected to be realizable in the laboratory within the next decade. 
We say "potential" because whether such simulations can be performed in a reasonable (from the perspective of 
the researcher or funding organization) amount of real time depends on how the implementation (for example, the 
vectorization, parallelization, memory usage) efficiently uses conventional computer resources. 

An important topic that is virtually unexplored is the implementation of fault-tolerant quantum computation 
on quantum computers other than the ideal one. The implementation of quantum error correction codes requires 
additional qubits. It is not known whether fault-tolerant quantum computation is robust with respect to systematic 
errors that inevitably affect the operation of physically realistic models of quantum computers. It most likely is 
not, and then it is good to know how to circumvent this problem. The complexity of this problem is such that it 
is not accessible to theoretical analysis (unless drastic, unrealistic simplifications are introduced), whereas computer 
simulation may give valuable insight. 

Decoherence is another topic that, albeit of much broader scope, is also very important for quantum computation. 
Although related to the issue of fault-tolerant quantum computation, the effect of the coupling between the quantum 
computer and its environment on the operation of the quantum computer is a problem that is hardly accessible to 
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FIG. 11: Microinstruction set window "grov-rf-k032" together with the quantum program windows "h_cor_l," "grovO," "grovl," 
and "grov2" that implement Grover's database search algorithm on a two-qubit NMR-like quantum computer for the cases 
go{x), gi{x), and g2{x) [see Eq. llUUt l. 



theoretical analysis, except perhaps for extremely idealized cases. Also here computer simulation may help. 

Finally, we believe that the software approach to quantum computation can help students gain insight into quantum 
physics during projects in which they improve their computational skills. The algorithms used to simulate quantum 
computers are fairly generic, that is, they solve differential equations of the parabolic type, and the concepts and 
techniques learned find applications in many other fields of computational science. 
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TABLE X: Programming languages for quantum computers. URLs last accessed on July 29, 2004. 



Name and 
latest release 



Description 



QCL (Beta) Version 0.5.1 
(0.6.1), March 30, 2004 



Q language Version 0.5.8, 
February 18, 2002 



QCL |1581 Il59l 116(3 . Il6l| (Quantum Computation Language) is a high-level, architecture- 
independent programming language for quantum computers, with a syntax derived from classical 
procedural languages like C or Pascal. Examples such as the quantum Fourier transform, Shor's 
algorithm, and Grover's algorithm are included. QCL has been developed under Linux; version 
0.5.0 compiles with the GNU C-f-|- compiler 3.3. The current version of QCL (sources and 1386 
Linux binaries) can be downloaded freely from: http://tph.tuwien.ac.at/~oemer/qcl.html 

The 



Q language jlOll . Il62l| is a C-|— I- implementation of a quantum programming language, 
source code can be downloaded freely from: http://sra.itc.it/people/serafini/qlang 



Quantum Superpositions 
Version 2.02, April 22, 2003 

QuBit, July 24, 2001 



Quantum Entanglement 
Version 0.32, June 5, 2002 



Q-gol Version 3, September 
11, 1998 



Quantum Fog Version 1.6, 
June 25, 2003 



QDD Version 0.2, February 
4, 2003 



Quantum Lambda Calculus 
2003 



Quantum Superpositions is a PERL library that enables programmers to use variables that 
can hold more than one value at the same time. The module can be downloaded freely from: 
http://search.cpan.org/~lembark/Quantum-Superpositions/lib/Quantum/Superpositions.pm 

QuBit is a C-l— I- library that supports Quantum Superpositions. The library is rewritten starting 
from the Quantum Superpositions library in PERL. A complete implementation of the QuBit 
code can be downloaded from: http://www.bluedust.com/qubit/default.asp 

Quantum Entanglement is a PERL library that allows the user to put variables in a superposition 
of states, have them in teract w ith each other, and then observe them. The module can be 
downloaded freely from: 'http://search.cpan.org/dist/Quantum-Entanglement 

Q-gol is a visual programming language. The fundamental symbols of the system are 
gates (represented by raised blocks), and directed wires. The graphical editor lets the 
user select gates, place them on sheets, and wire them together. An implementa- 
tion of Shor's algorithm is included. The source code can be downloaded freely from: 
|http://w ww. ifost.org.au/~gregb/q-gol/index. html 

Quantum Fog is a Macintosh application for modeling physical situations that exhibit 
quantum mechanical behavior. It is a tool for investigating and discussing quantum 
measurement problems graphically, in terms of quantum Bayesian nets. It simulates 
a general-purpose quantum computer. The software can be downloaded freely from: 
http://www.macupdate.com/info.php/id/12181 , http://www.ar-tiste.com 

QDD is a C-I--I- library that provides a relatively intuitive set of quantum computing constructs 
within the context of the C-I--I- programming environment. The emulation of quantum com- 
puting is based upon a Binary Decision Diagram representation of the quantum state. Shor's 
factoring algorithm is included in the QDD library, SHORNUF. The software can be downloaded 
freely from: http://thegreves.com/david/software 



Functional language based on Scheme for expressing and simulating quantum algorithms, 
at: http://www.het.brown.edu/people/andre/qlambda/index.htmls 



Access 



TABLE XI: Quantum compilers. URLs last accessed on July 29, 2004. 



Name and 
latest release 



Description 



Qubiter Version 1.1, April 6, 
1999 



GQC, 2002 



Qubiter is a quantum compiler written in C-I-+. Qubiter takes as input an arbitrary unitary 
matrix and returns as output an equivalent sequence of elementary quantum operations. The 
software can be downloaded freely from: http://www.ar-tiste.com/qubiter.html 



GQC I163II is an online quantum compiler that returns a circuit for the CNOT in terms of a 
user-specified unitary transformations. Access at: http://www.physics.uq.edu.au/gqc/ 



TABLE XII: Quantum circuit simulators. URLs last accessed on July 29, 2004. 



Name and 
latest release 



Description 



QCAD Version 1.80, May 8, 
2003 



QuaSi(2), March 18, 2002 



JaQuzzi Version 0.1, Jan- 
uary 14, 2001 



QCSimulator 

Version 1.1, March 8, 2000 



Libquantum Version 0.2.2, 
November 3, 2003 



OpenQUACS, May 22, 2000 



QuCalc Version 2.13, 
November 8, 2001 



QGAME Version 1, July 7, 
2002 



QCAD is a graphical Windows 98/2000 environment to design quantum circuits. It can export 
the circuits as BMP or EPS files, simulate the circuits, and show the states of the qubits. QCAD 
can be downloaded freely from: littp://acolyte.t. u-tokyo.ac.jp/~kaityo/qcad 

QuaSi(2) is a general-purpose quantum circuit simulator. It enables the user to build and 
simulate quantum circuits in a graphical user interface. Demo circuits for Shor's, Grover's and 
the Deutsch-Josza algorithm are included. QuaSi simulates up to 20 qubits. A Java applet is 
provided for use over the internet. The full version of QuaSi2 can be downloaded freely from: 
http://iaks-www.ira.uka.de/QIV/QuaSi/aboutquasi.html 



JaQuzzi 164IJ is an interactive quantum computer simulator to design, test, and visualize quan- 
tum algorithms with up to 20 qubits. The program can either run standalone or as a Web-based 
applet. To run jaQuzzi, a JAVA virtual machine of version 1.3 or higher is required. The 

software can be downloaded free ly from: 

'http://www.eng.buffalo.edu/~phygons/jaQuzzi/jaQuzzi.html' 

QCSimulator is a quantum computer simulator for Macintosh and Windows machines. A graph- 
ical user interface is used to build a circuit representation of a quantum algorithm and to 
simulate quantum algorithms by exchanging unitary elements with Mathematica Il65|. Shor's 
factorization algorithm, Grover's database search algorithm, the discrete Fourier transform, and 
an adder for two numbers are included. The complete software package can be ordered from: 
http://www.senko-corp.co.jp/qcs/ 

Libquantum is a C library for the simulation of an ideal quantum computer. Basic opera- 
tions for register manipulation such as the Hadamard gate or the ControUed-NOT gate are 
available. Measurements can be performed on either single qubits or the whole quantum regis- 
ter. Implementations of Shor's factoring algorithm, and Grover's search algorithm are included. 
Libquantum contains features to study decoherence and quantum error correction. Libquan- 
tum is developed on a GNU/Linux platform and requires the installation of a C compiler with 
complex number support. It can be downloaded freely from: [http://ww w.enyo.de/libquantumj 

OpenQUACS 'l66j (Open-Source QUAntum Computer Simulator) is a library written in Maple 
that simulates the capabilities of an ideal quantum computer. The simulator comes with a 
full tutorial. Several quantum algorithms such as Deutsch's algorithm, quantum teleportation, 
Grover's search algorithm and a quantum adder are included. The software can be downloaded 
freely from: http://userpages.umbc.edu/~cmccubl/quacs/quacs.html 



QCompute, July 1997 



Quantum, July 1997 



Eqcs Version 0.0.5, March 
19, 1999 



QCS, January 11, 2001 



QuCalc is a library of Mathematica functions [IBRl to simulate quantum circuits and solve 
problems of quantum computation. The Mathematica package can be downloaded freely from: 
http://crypto.cs.mcgill.ca/QuCalc 

QGAME |l67t Il68l | (Quantum Gate And Measurement Emulator) is a system, written in 
Common Lisp, that allows a user to run quantum computing algorithms on a digital com- 
puter. QGAME's graphical user interface (GUI) is a quick hack intended to allow peo- 
ple with no knowledge of Lisp to experiment with QGAME. It uses Macintosh Common 
Lisp (MCL) interface code and will work only under MacOS with MCL. Not all features 
of QGAME are available from the GUI. QGAME itself is platform independent (it will run 
on any platform for which a Common Lisp environment is available); only the GUI requires 
Macintosh Common Lisp. The full QGAME source code can be downloaded freely from: 
http://hampshire.edu/lspector/qgame.html 

QCompute is a quantum computer simulator written in Pascal that performs quantum gate 
operations on an arbitrary number of qubits. The source code and sample gate-networks for a 
1-bit and a 2-bit adder can be found in the appendices to a thesis that can be downloaded from: 
http://w3.physics.uiuc.edu/~menscher/quantum.html 

Quantum is a quantum circuit simulator written in C-|— f for a Windows environment. It in- 
cludes an implementation of Shor's algorithm. The software can be downloaded freely from: 
http://www.themilkyway.com/quantum 

Eqcs is a library allowing clients to simulate a quantum computer. It includes a pro- 
gram showing the creation of a CNOT gate. The software can be downloaded freely from: 
http://home.snafu.de/pbelkner/eqcs/index.html 

QCS (Quantum Computer Simulator) is a quantum computer library in C-I--I-. The software can 
be downloaded freely from: http://www-imai. is. s. u-tokyo.ac.jp/~tokunaga/QCS/simulator.htmll 



TABLE XIII: Quantum Computer Emulators. URLs last accessed on July 29, 2004. 



Name and 
latest release 



Description 



QCE Version 8.1.1, June 27, 
2004 



QSS, June 2000 



QCE 123, llSol (Quantum Computer Emulator) is a software tool that emulates ideal quantum 
computers as well as physical implementations of quantum computer hardware. QCE uses time- 
dependent Hamiltonians and unitary time evolutions to simulate the physical processes that 
govern the operation of a hardware quantum processor. QCE provides an environment to debug 
and execute quantum algorithms under realistic experimental conditions. The QCE package 
includes ma ny ex amples such as Shor's algorithm for factoring integers, order finding, number 
partitioning |102| . the quantum Fourier transform, various implementations of the Deutsch-Josza 
algorithm, and Grover's database search on ideal and more realistic quantum computers, such 
as those used in the 2-qubit NMR quantum computer. The software consists of a Graphical 
User Interface and the simulator itself. It runs under Windows 98/NT4/2000/ME/(XP with 
SPl) and Linux+VMware on Intel/ AMD processor machines. The software can be downloaded 
freely from: http://www.compphys.org/qce.htm 

QSS 169] (Quantum System Simulator) is a software tool that simulates quantum computations 
with time-dependent Hamiltonians. It provides a simple and convenient graphical user interface 
that enables users to specify complex Hamiltonians as sums of smaller ones. The simulator 
engine is designed as a self-contained module, so that it is independent of the user interface, and 
can be easily enhanced. The QSS runs on a Windows operating system and can be downloaded 
freely from: http://web.mit.edu/scottsch/www/qc 



TABLE XIV: Pedagogical software. URLs last accessed on July 29, 2004. 



Name and 
latest release 



Description 



Quantum Turing machine 
simulator 2002 



QTM simulator, 1995 



Quantum Search Simulator, 
October 10, 2002 



Grover's algorithm, 
June 2001 

CS 596 Quantum Comput- 
ing, Spring 1999 



Shor's algorithm, June 2001 



Shor's algorithm simulator, 
January 23, 2002 



Mathematica toolkit to construct, run, and research quantum Turing ma- 
chin es. S ubscribers to The Mathematica Journal can download the toolkit from: 

http://www.mathematica-journal.com/issue/v8i3/features/hertel 

Classical simulator of a quantum Turing machine written in C++. The machine has one one- 
sided infinite tape and one read/write-head and simulates a Fourier transform on a parity 
function. It can be downloaded freely from: http://www.lri.fr/~durr/Attic/qtm 

Quantum Search Simulator is a Java applet that demonstrates the operation of Grover's quan- 
tum search algorithm on a database of four items on a quantum computer based on optical 
interference. Access at: http://strc.herts.ac.uk/tp/info/qucomp/qucompApplet.html 

Mathematica-compatible notebook demonstrating Grover's algorithm. It can be downloaded 
freely from: http://www.cs.caltech.edu/"' chenyang/cs20proj-grover6.nb 

Matlab programs that demonstrate some features of quantum computing. Demonstrations of 
the RSA algorithm and Shor's algorithm on a conventional computer are included. It can be 
downloaded freely from: http://www.sci.sdsu.edu/Faculty/Don.Short/QuantumC/cs662.htm 

QCL code demonstrating Shor's algorithm. It can be downloaded freely from: 
http://www.cs.caltech.edu/"' chenyang/myshor.qcl 

C-I--I- program that simulates the operation of a quantum computer performing Shor's algorithm. 
It can be downloaded freely from: http://alumni.imsa.edu/~matth/quant 
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